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CHAPTER  1 


INTRODnCTION 


The  chained  aggregation  procedure  was  first  introduced  in  [1]  in  the 
context  of  reduced— order  modelling  of  large-scale  systems.  The  chained 
aggregation  procedure  transforms  the  original  system  into  the  Generalized 
Hessenberg  Representation  (GHR).  Once  placed  in  the  GHR,  the  reduced-order 
modelling  analysis  is  carried  out. 

Since  the  introduction  of  the  chained  aggregation  procedure,  much  more 
research  has  been  done  [2,3,4,51.  The  work  that  has  followed  has  not  been 
constrained  to  the  reduced-order  modelling  problem.  But,  it  has  expanded 
the  possible  applications  of  the  chained  aggregation  procedure  to  include 
many  control  system  problems.  For  a  thorough  discussion  of  chained 
aggregation  and  the  (SR  for  control  system  design,  the  author  recommends 
consulting  [6],  where  this  issue  has  been  specifically  addressed  in  a 
geometric  setting. 

The  algorithm  developed  in  this  thesis  uses  the  numerical  advantages 
associated  with  orthogonal  matrices  to  implement  the  basic  chained 
aggregation  procedure.  Included  in  the  algorithm  is  an  enhanced  procedure, 
called  modified  chained  aggregation  (MCA),  which  has  been  introduced  in 
earlier  literature  [2,6,7].  Several  numerical  advantages  of  orthogonal 
autrices  are  included  herein;  more  discussion  may  be  found  in 


[8,9,10,11,12] . 
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The  orthogonml  matrices  are  obtained  during  the  algorithm  by  using  the 
singular  value  decomposition  (SVD)  of  particular  matrices.  Briefly,  the  SVD 
is  used  to  identify  the  linearly  independent  and  linearly  dependent  rows  and 
columns  of  the  particular  matrix.  With  the  identification,  the 
decomposition  generates  orthonormal  vectors  spanning  the  subspaces 
associated  with  the  above  sets  of  rows  and  columns.  These  orthonormal 
vectors  may  then  be  judiciously  grouped  to  form  orthogonal  matrices  to  be 
used  within  the  algorithm. 

The  organization  of  this  thesis  is  as  follows.  Chapter  2  presents  the 
SVD  and  steps  through  the  chained  aggregation  algorithm  showing  how  SVD  has 
been  incorporated.  The  development  of  the  modified  chained  aggregation 
algorithm  is  contained  in  Chapter  3.  The  chained  aggregation  algorithm  is 
used  in  Chapter  4  to  help  in  the  design  of  an  output  feedback  controller  for 
a  particular  large  flexible  space  structure.  Chapter  S  presents  the 
conclusion. 

A  summary  of  definitions  and  the  notation  used  throughout  this  thesis 
follow.  Given  a  subspace  S,  denoted  by  bold,  capital  roman  letters,  its 
orthogonal  complement  will  be  represented  by  S  .  The  range  space  and  the 
null  space  of  a  matrix  will  be  denoted  by  S[>]  and  N[*],  respectively. 
is  used  to  represent  the  pseudo-inverse  of  the  matrix  A.  The  transpose  of  a 
matrix  A  will  be  denoted  by  A  .  The  set  of  all  mxn  matrices  of  rank  r  with 
coefficients  in  the  real  number  field.  R.  will  be  denoted  by  The 

working  precision  of  the  computer  will  be  represented  by  e.  On  a  given 
computer,  the  value  of  e  equals  the  smallest  number  which  when  added  to  one 
equals  one.  The  singular  values  of  a  matrix  will  be  denoted  by  o^.  The 
norms  I  I *112  and  I  I "lip  correspond  to  the  matrix  2-norm  and  the 


matrix  Frobenins  norm,  respectively.  If  A  e  R®***,  then 


llAll^ 


maximize  ,  ^  , 

ly^AxI  . 

IUIl2=llyll2  =  1 


where 


ll^llj  =  (  Z  =  (xTx)l/2 

i^l 

ffi  n 

(  I  Z  af.)l/2  . 

i=l  j=l 

The  dimension  of  a  vector  space  will  be  denoted  by  d(«).  The  notation 


X  = 

llAlIp  = 


represents  the  vector  space  spanned  by  the  non-zero  elements  of  the  columns 
of  X  as  they  vary  over  the  real  numbers.  The  rank  of  a  matrix  will  be 
denoted  by  p[ • ] . 
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CHAPTER  2 


SINGULAR  VALUE  DECOBJPOSITION  AND  THE  CHAINED  AGGREGATION  ALGORITHM 


2.1,.  Introduction 

Chained  aggregation  as  introduced  in  [1]  identifies  the  information 
structure  of  the  system  by  aggregating  the  system  with  respect  to  the 
output.  Once  aggregated,  the  system  exhibits  the  Generalized  Hessenberg 
Representation  (GHR)  structure  which  has  been  shown  to  st'  ..  te  many 
possible  design  procedures  [1,2,3 ,4,5,6] . 

A  design  procedure  used  in  this  thesis  is  model  rednct  or  Model 
reduction  arises  out  of  the  simple  interconnecting  structure  displayed  in 
the  GHR  and  by  identifying  the  feedback  coupling  between  subsystems  which  is 
weak  or  nonexistent.  The  resulting  model  retains  the  strongly  observable 
modes  from  the  outputs  rather  than  retaining  the  dominant  system  modes  as  in 
modal  reduction.  Throughout  the  following  discussion  the  subsystem  composed 
of  the  strongly  observable  modes  will  be  referred  to  as  the  aggregate 
subsystem,  while  the  residual  subsystem  will  refer  to  the  remaining 
subsystem. 

Associated  with  the  chained  aggregation  procedure  are  many  state  space 
transformations.  It  is  shown  below  these  transformations  can  be  performed 
using  orthogonal  matrices  which  are  numerically  robust,  in  contrast  to  the 
nonorthogonal  transformations  described  in  most  of  the  previous  literature 
[1,2, 3, 4, 7] . 


r 


The  following  section  reviews  the  singnl°' 


iiiue  decomposition  ( SVD) 


along  with  some  of  'ts  properties.  The  SVD  is  then  incorporated  into  the 
steps  of  the  chained  aggregation  algorithm  as  shown  in  Section  3.  Section  4 
highlights  a  consequence  of  the  chained  aggregation  algorithm,  its  ability 
to  check  the  observability  and  controllability  of  a  system. 

2.2.  Singular  Valne  Decomposition  and  Its  Propert ies 

The  SVD  will  be  introduced  by  the  following  theorem. 

Theorem  2.2.1[8]  :  Let  A  €  Then  there  exist  orthogonal  matrices 

U  e  and  V  €  R'^®  such  that 


A  =  U  Z  v'^ 


(2.2.1) 


where 


I 


S  0 
0  0 


and  S  -  diBg(<j^,<,2, . . .  .Of)  with  .  .^a^yO. 

Proof :  See  [8] . 

T 

The  product  U  Z  V  is  the  singular  value  decomposition  of  the  A  matrix. 

The  numbers  ,02 , . . .  ,<»y  together  with  are  called  the  singular 

T 

values  of  A  and  are  the  positive  square  roots  of  the  eigenvalues  of  A  A. 

The  following  well-known  properties  of  SVD  make  it  useful  in  the 
chained  aggregation  algorithm. 
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Two  different  matrix  norms  associated  with  can  easily  be 

defined: 


llAlIj  =  ai 

llAlIp  =  +  ...  +  . 

Our  main  interest  in  singular  values  will  be  for  rank  determination. 
Thus,  knowing  that  the  singular  values  are  not  very  sensitive  to 
perturbations  in  the  matrix  insures  a  good  rank  determination.  The  change  in 
the  singular  values  are  known  to  be  bounded  by  the  magnitude  of  the  matrix 
perturbation. 

Theorem  2.2  ;  Let  A,B  6  R®**'  have  singular  values  ® 1  . . .  >.  and 

2.  2.  respectively.  Then 

-  T^l  <  II  A  -  B  II2  (i  =  l,2,...,n). 


Proof:  See  [11] . 

The  concern  over  matrix  perturbations  is  related  to  the  fact  that  infinite 
precision  arithmetic  cannot  be  performed  on  a  computer.  For  this  reason, 
when  the  singular  values  of  a  matrix  A  are  desired,  the  computed  singular 
values  are  actually  the  exact  singular  values  of  a  matrix  slightly  perturbed 
from  A,  say  A  E.  A  more  extensive  discussion  can  be  found  in  [8,10] . 

To  determine  the  rank  of  a  matrix,  all  of  its  singular  values  are 
compared  with  y  =  e«||A||,  where  the  particular  norm  used  can  be  selected  by 


7 


the  user;  the  number  of  ®^'s  greater  than  y  detemines  the  rank  [10]. 
Looked  at  in  this  way.  the  smallest  singnlar  value,  a^,  greater  that  y  gives 
a  measure  as  to  how  far  away  the  A  matrix  is  from  another  matrix  of  smaller 
rank. 


Another  important  outcome  of  SVD  is  the  identification  of  the  fonr 
fundamental  subspaces  [9]  of  a  matrix.  Suppose  F  €  has  an  SVD  given  by 


F 


D  I 


- 

s  o' 

’vl’ 

°1  '  «2 

0  0 

Vi 

(2.2.2) 


where  S  =  diag(<T^, , . ,  ^0^)  i^ith  Oj  >_  ...  >.  Oj  >  0  and  U  and  V  are  partitioned 

compatably,  i.e.,  Uj  e  e  etc.  With  this  notation  Dj,  C2  and 

Vf.  V2  produce  orthonormal  bases  for  the  four  fundamental  subspaces,  £(F] , 

X  X 

£  [F] ,  N  [F] ,  N[F]  [8].  Figure  2.1  redrawn  from  [8]  relates  these 

subspaces . 

Associated  with  the  matrix  F  €  R*****^  are  two  different  vector  spaces,  R*" 

and  R°.  Figure  2.1  illustrates  how  the  F  matrix  can  induce  direct  sum 

decompositions  of  these  two  vector  spaces.  In  R°  the  columns  of  V2  form  an 
orthogonal  basis  for  the  N[F]  while  the  columns  of  span  the  remainder  of 
R*^  with  an  orthogonal  basis.  Similarly,  the  columns  of  form  an 

orthogonal  basis  for  the  £[F]  in  R™  while  the  columns  of  D2  span  the 
remainder  of  R"*  with  an  orthogonal  basis.  The  mappings  between  these 

subspaces  is  also  indicated.  The  F  matrix  maps  the  space  spanned  by  into 

the  £[F] ,  i.e.,  the  space  spanned  by  Dj^.  The  pseudo-inverse  of  F  performs  a 

map  in  the  opposite  direction.  The  space  spanned  by  the  columns  of  V2  is 
mapped  into  the  origin  in  R*"  by  the  F  matrix.  The  space  spanned  by  the 
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colunns  of  is  mapped  into  the  origin  in  R®  by  F^. 

The  pseudo-inverse  of  F  can  be  obtained  from  the  SVD  of  F 


V  i‘*‘ 


(2.2.3) 


where  S“^  =  diag(i, . . .  >  . . .  2  >  0,  and  D  and  V  are  obtained 

from  (2.2.2)  . 

In  general,  SVD  is  the  only  numerically  reliable  method  of  generating 
basis  for  these  subspaces,  since  U  and  V  are  partitioned  into  [U^  ,  |}2]  and 
,  V2] ,  according  to  the  smallest  singular  value  [8]. 

Further  discussion  of  Figure  2.1  will  be  postponed  until  after  the 
chained  aggregation  algorithm  has  been  presented. 

The  advantage  of  using  SVD  within  the  chained  aggregation  algorithm 
lies  in  the  fact  that  the  bases  vectors  generated  by  SVD  are  orthonormal. 
Thus,  constructing  transformation  matrices  from  these  orthonormal  vectors 
results  in  orthogonal  transformations.  Several  well-known  numerical 
advantages  associated  with  orthogonal  matrices  useful  in  the  chained 
aggregation  algorithm  are: 

(1)  Orthogonal  matrices  are  easy  to  invert,  IT^  = 

(2)  Orthogonal  matrices  are  perfectly  conditioned  with  respect  to  the 
2-norm,  IIADII^  =  IIAII2. 


Ji 


(3)  Orthogonal  matrices  lend  themselves  to  backward  error  analyses 
[11].  For  example,  suppose  an  error  F  is  introduced  into  the  result  of  an 
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orthogoaal  trsnsforaation.  Let  E  =  DFD^.  Fron  the  second  characteristic 
above 


MEII2  =  llDF0^ll2  =  l|FII^Il2  =  IIFII2 


and 


U^(A  +  E)D  =  O'^AO  +  D^D  =  D^AD  +  F  . 

In  other  words,  a  perturbation  in  the  result  can  be  accounted  for  by  a 
perturbation  of  the  same  magnitude  in  the  original  problem.  This  guarantees 
that  if  there  exists  an  uncertainty  in  the  original  A  matrix  of  magnitude  6, 
then  the  resulting  transformed  system  will  have  an  uncertainty  of  the  same  & 
magnitude. 

(4)  Orthogonal  transformations  only  rotate  the  axes,  each  axis 
maintaining  its  exact  relation  to  the  others  throughout  [12j .  It  should  be 
mentioned  that  not  all  orthogonal  matrices  represent  pure  rotations  as  can 
be  seen  from  a  simple  example.  Consider  the  orthogonal  matrix 


This  matrix  reflects  every  point  (x,y)  into  its  mirror  image  (y.x)  across 
0 

the  45  line  y  =  x,  which  is  not  a  true  rotation. 

In  [1]  orthogonal  matrices  have  not  been  used  for  the  transformation 


matrices  and  consequently  the  matrices  used  can  be  very  ill  conditioned  and 
result  in  a  numerically  unstable  algorithm. 
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Th*  Chaingd  Aagreaation  Alaoritha 

The  chained  aggregation  algorithm  transforms  any  given  system  of  the 

form 


X  >•  Ax  +  Bu 


(2.3.1a) 


y  =  Cx  ,  (2.3.1b) 


rixn 

where  A  €  B  €  and  C  6  R  into  the  6HR  nsing  only  the 

information  strnctnre  of  the  system.  The  GHR  is  given  by 


e 

z  Fz  +  Gu  (2.3.2a) 


y  =  Dz  ,  (2.3.2b) 


where 


J 
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steps  comprise  the  chslned  s(gregatioo  algorithm. 


Alaoritl 


r-ixn 

Let  A  e  B  e  R““,  and  C  6  R  ^  with  C  having  full  row  rank. 
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Step  Perform  a  colnmn  compression  (101  on  the  C  matrix  by  nsing  SVD.  The 
linearly  independent  columns  of  C  will  be  compressed  to  the  left  by  post 
mnltlplying  C  by  the  V  matrix  resulting  from  SVD.  If  V*  *  [V^  ,  v|]  with 

nxr-i  nx(n~ri) 

€  R  ^  and  V|  €  R  ^  ,  then 

c  = 


resnlts  in 

CV^  =  C[VJ  ,  V^l  =  (Dj  0]  . 
Step  1".  Initialize  the  transformation  matrix 


-j-i  _  yiT  ^ 

Step  3.:  Perform  a  state  space  transformation  on  the  matrix.  This 
transformation  can  be  broken  down  into  four  smaller  and  separate 
transformations.  Since  at  most  only  two  of  the  four  smaller  transformations 
are  needed,  this  allows  for  fewer  calculations  during  the  execution  of  the 
algorithm. 


iTAiyi 

11 


iT.ivi 

2  1 


Vf  Aiyi  - 
2  2  -i 


Afl  Ai2 

21  ^22 


Let 
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Compute  only  ®  ®  ^  *tep. 


Step  4:  Compress  the  columns  of  A^2  °sing  SVD. 


Ai^  =  ui+lj;i+lv(i+l)T 


^^i,i+l  ' 


where  F.  G  and 


Step  4:  Check  the  singular  values  of  A22* 

If  all  ®j's  <  Y  (j  =  l,2,...,r^),  then  exit  1. 

If  all  a^.'s  >  y  (j  =  l,2,...,r|^2  ^  r^) ,  then  exit  2. 
Otherwise,  continue  with  r  -  t  +  (number  of  Oj's  >.  y  ) 


Step  6.:  Update  the  transformation  matrix. 


•j'i+l 


Q  yC  i+1)  T 


where  implies  I  €  R 


TXT 


Step  2‘  Calculate  the  A22  submatrix. 


"22 


v*'^A‘v‘  e  rp*p 


Step  8:  Let  A^"*"^  =  A^. ,  i  =  i  +  1.  Go  to  Step  3. 


22 


Exit  1:  The  system  aggregates.  F  =  TAT^,  G  =  TB,  D  =  CT^,  where  F^  =  0. 


Exit  2:  The  system  does  not  aggregate.  F  =  TAT^,  G  =  TB,  D  =  CT^,  where 


J 


IS 


full  column  rank. 


Upon  exit  from  tkis  algorithm  the  original  system  has  been  transformed  into 
the  GES  nsing  orthogonal  matrices,  a  numerically  stable  transformation 
process . 


The  chained  aggregation  algorithm  identifies  the  supremal  A-invariant 
subspace  in  the  N[C]  [13] .  This  can  be  understood  by  referring  to  Eq. 
(2.3.2).  In  the  new  bases  the  supremal  A-invariant  snbspace  in  the  N[C]  is 
immediately  seen  to  be 


sp 


where  X  €  R 


szl 


and  s  =  d(F^  . 


Any  vector  lying  in  this  space  is  obviously  A-invariant,  i.e.. 


and  lies  in  the  N[C] , 


IDj  0  0] 


=  [0]. 


To  better  understand  Figure  2.1,  consider  the  system  at  step  3  on  the 
first  pass  through  the  algorithm.  Aggregation  will  occur  at  this  stage  if 
and  only  if  the  N[C]  is  A-invariant,  i.e.,  AN[C]C  N[C] .  If  the  F  matrix  in 
Figure  2.1  is  replaced  by  the  C  matrix  of  the  system,  the  columns  of  are 
seen  to  span  the  M[C]  and  therefore  AV^C  V2.  This  immediately  forces  the 


snbmatrix  to  be  zero 
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Ai2  =  vJaV^C  vJVj  =  0  . 

Thus,  by  using  Figure  2.1  and  noting  what  causes  aggregation,  an  understanding 
of  why  the  submatrix  must  be  zero  can  be  obtained. 

^.4.  Observability  and  Controllability  Tests 

An  immediate  consequence  of  the  algorithm  is  its  ability  to  test  the 

observability  of  a  system  in  a  numerically  stable  manner.  By  considering  the 
T  T 

dual  problem  (A  ,  B  ).  the  controllability  of  a  system  can  also  be  checked. 
Algorithms  similar  to  the  chained  aggregation  algorithm  dealing  specifically 
with  the  controllability  problem  have  been  given  in  [10,12,14]. 

The  following  definition  is  needed. 

Definition  2.£.1:  The  pair  (A,C)  is  said  to  aggregate,  if  when  represented 
in  the  CHR  the  j.  snbmatriz  equals  zero. 

The  tests  using  the  chained  aggregation  algorithm  may  now  be  given. 

Theorem  2.4.2:  The  following  statements  are  equivalent: 

1.  The  pair  (A,C)  is  unobservable. 

2.  The  system  (A,C)  will  aggregate. 

Proof :  (  1  implies  2  ). 

If  the  pair  (A,C)  is  unobservable,  then  the  Hautus  test  will  result  in 
a  matrix  of  rank  less  than  n.  Since  observability  is  invariant  to  a  state 
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space  transformation,  the  Bantus  test  may  be  applied  to  the  GBR  induced  by 
the  pair  (A,C) 


Because  of  the  GBR  construction,  the  only  manner  in  which  the  above  matrix 
can  have  rank  less  than  n  for  any  X  is  for  the  ^  submatrix  to  be  zero. 
This  is  exactly  the  condition  for  aggregation. 

(  2  implies  1  )  is  now  clear.  If  the  system  aggregates,  then  F  .  v  ~  ^ 
and  the  pair  (A,C)  is  unobservable,  since  the  rank  of  the  above  matrix  is 
less  than  n  for  any  eigenvalue  of  the  system.  □ 

The  dual  result  follows  immediately. 

Theorem  2.4.3:  The  following  statements  are  equivalent: 


18 


1.  The  pair  (A,B)  is  ancoatrol lable . 

2.  The  system  (A^  .  B^)  will  aggregate. 

Proof :  Similar  to  the  above  proof. 


CHAPTER  3 


MODIFIED  CHAINED  AGGREGATION 

3..1 .  Introdoct  ion 

In  the  basic  chained  aggregation  algorithm,  only  the  output  information 
structure  of  the  system  is  used.  By  incorporating  the  input  structure,  the 
system  can  be  forced  to  aggregate  ,  using  state  feedback,  when  the  R[B^] 
contains  the  RIAj^^J  ,  where  Aj2  snd  Bj  are  the  submatrices  generated  during 
one  of  the  steps  in  the  algorithm.  The  ability  of  the  input  to  satisfy  the 
above  condition  is  determined  using  the  singular  value  decomposition  (SVD) 
to  identify  the  rows  of  Bj^  which  are  linearly  independent.  By  the  proper 
choice  of  input,  u,  these  linearly  independent  rows  can  then  be  used  to 
annihilate  the  largest  possible  subblock  of  the  A^^^  submatrix.  Further 
explanation  is  given  below. 

This  process,  called  modified  chained  aggregation  (MCA),  has  been 
developed  in  [2]  for  use  in  Three-Control-Component-Design  (TCCD) .  The 
concept  of  TCCD  [2,6,7]  arises  in  the  Generalized  Hessenberg  Representation 
(GHR)  structure.  Briefly,  this  hierarchical  design  procedure  uses  the  input 
in  three  specific  ways.  First,  the  input  is  used  to  force  aggregation. 
Second,  the  aggregate  dynamics  are  adjusted  to  the  specifications  of  the 
designer.  Third,  if  enough  input  structure  exists  the  residual  dynamics  can 
be  adjusted. 

The  application  of  the  above  concepts  has  been  carried  over  into 
controller  design  for  a  special  class  of  nonlinear  systems  [S],  thus 


J 
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broadening  the  scope  of  problems  capable  of  solution  using  the  developed 
algorithms . 

An  algorithm  [10]  similar  to  MCA  appearing  at  approximately  the  same 
time,  although  the  implementation  had  not  been  completed,  was  motivated  by 
determination  of  the  supremal  ^A,B) -invariant  snbspace  in  the  N[C]  [13]. 

Section  2  demonstrates  how  the  input  structure  is  identified  and  used 
to  enhance  the  chained  aggregation  algorithm.  The  steps  of  an  implemented 
algorithm  which  accomplishes  the  goals  of  both  procedures  above  [2,10]  are 
contained  in  Section  3. 

^.2.  Identifying  The  Input  Structure 

To  motivate  the  use  of  the  input  structure  within  the  chained 
aggregation  algorithm,  consider  the  following  example. 

Suppose  the  system  after  one  stage  of  chained  aggregation  has  the 
following  form: 


”  *“ 

—  “ 

*1 

All 

Ai2 

Bl 

^21 

^22 

*2 
_  * 

«2 

(3.2.1a) 


y 


(3.2.1b) 


Aggregation  occurs  when  A^^  =  o.  If  p  0.  but  R[Ai2K  l[Bj]  ,  then  by- 
selecting  the  input  u  =  -  system  can  be  forced  to  aggregate. 
This  assumes  the  states  are  available  for  feedback.  If  these  states  are 


21 


not  explicitly  available,  then  some  t3rpe  of  dynamic  feedback  mast  be 
introduced  to  reconstruct  them.  If  ,  then  the  MCA  algorithm 

identifies  the  largest  subspace  of  A^^  contained  in  the  and  performs 

another  step  of  chained  aggregation  using  the  subspace  of  A^^  not  contained 
in  the  E[Bj]  as  the  aggregating  matrix. 

SVD  is  used  to  calculate  the  K[B^] ,  thus  a  numerically  stable 
computation  is  obtained.  To  identify  the  EIB^^]  a  row  compression  is 
performed  within  the  MCA  algorithm.  The  orthogonal  transformation  which 
achieves  this  row  compression  is  then  treated  as  a  state  space 

transformation  which  is  used  to  transform  the  system  matrix.  After 

transforming  the  system  matrix,  the  snbmatrix  is  divided  into  two 

submatrices;  one  submatrix  A^^  is  not  in  the  RtB^] ,  and  the  other  snbmatrix 

If  ^1  represents  the  compressed  rows  of  Bj,  then  by  using  the 

feedback  u  =  -  f^e  k-^2  submatrix  can  be  annihilated  in  the  system 

matrix.  Again,  the  states  have  been  assumed  to  be  available. 


Consider  the  SVD  of  the  Bj^  submatrix 


U  I  = 


.  0^ 


S  0 
0  0 


(3.2.2) 


n.  X  r  n, 

where  D^  e  R  .  G  R 

are  the  singular  values  of 

rows  of  the  B^^  snbmatrix  are 

must  be  premultiplied  by  [U^ 


X  (uj-t) 

,  S  =  diag(a2,a2> . . .  .Oj.)  .  and  the  ® j  '  s 
For  ease  in  implementing  the  algorithm  the 
compressed  down.  To  achieve  this  result,  B^ 
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R1 

B,  = 

1 

a 

CO 

_ 1 

r  T 1 

yl. 

1 

flh. 

yl. 

h 


l)[Uj 


-  — 

T 

svi 

"I”2 

0 

(3.2.3) 


The  resulting  submatrix  has  full  row  rank.  If  the  original  B  matrix  had 
the  form  B  =  [B^  ,  B^]^,  then  the  desired  row  compression  in  the  B^^ 
submatrix  can  be  carried  out  by  performing  a  state  space  transformation  with 
the  following  matrix. 


T 


(3.2.4) 


Performing  a  state  space  transformation  with  this  matrix  results  in  the 
following  structure: 


tat’^ 


1 — 
O 

r  I 

T 

^11  Aj2 

U2  0 

0 

-  ®  K 

^21  ^22 
_ 

- 1 

1-4 

o 

o 

_ J 
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U2^11% 

02^11% 

02^12 

A'll 

Ai2 

«1^11°2 

Ul^l°l 

Di^12 

= 

A ' 

Ai2 

_  ■'^12^2 

^21*^1 

^22  _ 

21 

- 1 

< 

(3. 2., 5) 


If  ^12  ^  If  as«d  as  the  aggregating  matrix  in  the  next  step  of 

chained  aggregation. 


The  above  technique  is  nsed  throughout  the  MCA  algorithm. 


3,.^.  The  Modified  Chained  Aggregation  A1  finrithm 


,  ,  r.xn 

Initialization;  =  A  €  =  3  €  R“®.  C  €  R  . 

p=n,  T=0,  0=0,  i=l. 


=  0,  q  =  0, 


Step  _1:  Compress  the  columns  of  C  using  SVD 


C  =  0^ 


This  implies 


yi  yi 
1  '  '^2 


'  ^2 


0 


0  0 


D][S‘  ,  0 


nzri 


where  €  R  ^  and  €  R 


nx(n-rn ) 


.  Let  a  r. 


Step  2:  Initialize  the  transformation  matrix 


•j*  I  =  y 


Step  Calculate  the  A^^  submatrix 


^12  =  vj^AivI  ,  aJj  €  R«(p-“) 
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If  A^2  ~  0,  then  exit  1. 


Step  4:  Cxlcnlete 


A,  10. 

^  .  B^ 

0 


TiT»  ''i  Pl“  A.  (n-p,-n)xai  * 

where  0  6  R’'“.  B^  e  R  .  end  B^  €  R  with  =  o  +  p 


A  1 

Step  Conpress  the  rows  of  Bj  down  using  SVD 


Bj  =  fii  ii  viT  . 


of 

A 

fixT 


^  wT^  i 

oJ^b; 


where  €  R^™  . 

Note:  The  veloe  for  B  has  changed  between  steps  4  and  S.  It  now 

Ai 

equals  the  nnaber  of  linearly  independent  rows  in  the  B^  snbmatrix. 


Step  b:  Update  the  B  natrix 


2S 


.i+1 


0 

0 


L®2J 


(Tizm) 

((Pj-P)xin) 

(pzin) 


( {n“^ j— n) xm) 

Step  2=  Update  the  transformation  matrix 


<1  _ 


02 

^iX 

Dl 

0 


0 

0 

0 

[ 

p-a 


Ti 


Step  Update  the  submatrix 


*12  '  »f*i2 


,  where  Aj2  ^  B 


(pj-p)x(p-a) 


3^22  “  ^12  ®  . 


If  A^2  *  0,  then  exit  1. 


Step  £:  Calculate  the  A22  submatrix 


AI2  =  viTAiyi  .  where  A^,  €  r( P-<i) x( p-a) 


'2  2 
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p-a 


■z  =  T  +  a 


1  = 


Step  10 :  Compress  the  columns  of 


i  +  1 
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Ai-1  =  oi  r*  viT  . 


^i-1 

^12 


ui 

Vi  , 


D‘S^  .  0 


where  6  £***“  and  6  rP^^P"®) 


If  a  »  p.  then  exit  2, 


Step  11;  Update  the  transformation  matrix 


It  0 
0 

0  viT 


ri““l 


Step  12:  Compute  the  new  submatrix 


7(  i-1)  vi 
^12  ^^2 

viT.i  yi 
1  *22^2 


where 


If  Aj2  =  0»  then  exit  1. 


Step  13 :  n  »  t  -  a  .  Go  to  Step  4, 


Exit  1,:  TTie  system  will  aggregate.  F  =  TAt'^,  G  =  TB,  and  D  =  CT^,  where 
T  =  T^. 

Exit  2:  The  system  will  not  aggregate.  F  =  TAT^,  G  =  TB,  and  D  =  Ct"^,  where 
T  =  T‘. 
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CHAPTER  4 


OUTPUT  FEEDBACK  DESIGN  USING  THE  CHAINED  AGGREGATION  ALGORITHM 


4.1.  Introdact ion 

Becnnse  of  the  conpotational  and  numerical  difficnlties  associated  with 
large  scale  systems,  a  reduced-order  model  of  the  system  is  advantageous 
during  the  design  process.  As  stated  in  Chapter  2,  the  Generalized 
Hessenberg  Representation  (GHR)  of  a  system  lends  itself  nicely  to  the 
possibility  of  identifying  a  suitable  reduced-order  model.  By  comparing  the 
sizes  of  the  submatrices  using  an  appropriately  selected  norm  generated 
during  each  stage  of  chained  aggregation,  a  trade-off  between  reduced-order 
model  dimension  and  subsystem  coupling  can  be  made  to  obtain  the  desired 
reduced-order  model.  The  resulting  reduced-order  model  contains  the 
strongly  observable  modes,  by  construction,  which  may  or  may  not  be  the 
dominant  system  modes. 

This  technique  of  generating  a  reduced-order  model  has  been  conducted 
on  a  particular  large  space  structure.  The  general  space  structure  problem 
is  introduced  in  Section  2.  Section  3  contains  the  description  of  the 
particular  structure  which  was  studied.  The  design  procedure  which  was 
carried  out  is  detailed  in  the  final  section. 

4.2.  Introduction  of  the  General  Larae  Space  Structure  Problem 

Large  Space  Structure  (LSS)  problems  have  received  a  great  deal  of 
attention.  The  LSS  themselves  are  usually  quite  flexible  because  of  weight 
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restrictions  isiposed  during  their  transport  or  deploTiaent  in  space.  The 
problens  are  farther  coaiplicated  due  to  the  zero  damping  environment  of 
space  and  the  light  natural  damping  of  the  structure  itself.  A  recent 
review  of  the  literature  on  LSS  Control  can  be  found  in  [IS]  with  more 
information  available  in  [16-23]  . 

Initially  these  systems  are  often  modelled  by  partial  differential 
equations.  For  a  practical  solution,  one  must  reduce  the  infinite 
dimensional  problem  down  to  one  of  finite  dimension.  The  most  popular 
method  of  reducing  the  infinite  dimensional  problem  has  been  the  finite 
element  method  [15],  a  structural  analysis  technique.  The  differential 
equations  resulting  from  the  finite  element  method  are 


Mq  +  Kq  *  Bpu 


(4.2.18) 


y  =  C^q  ,  (4.2.1b) 

where  q  €  R“,  u  €  R®,  y  €  R^,  and  the  constant  matrices  M,  K,  Bp,  and  are 
of  cor.patable  dimensions.  N  is  the  mass  matrix,  which  in  general  is 
positive  semidef inite ,  and  K  is  the  stiffness  matrix,  which  is  positive 
definite.  The  physical  nature  of  most  LSS  problems  allows  the  damping  to  be 
neglected  in  the  modelling,  i.e.,  there  are  no  terms  involving  q  in  Eq. 
(4.2.1a) . 

To  perform  chained  aggregation  on  the  system  it  must  be  represented  in 


state  space  form.  A  common  method  to  obtain  this  representation  has  been  to 
simultaneously  diagonalize  M  and  K  with  a  unitary  matrix  i  [17,18]  such  that 
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«^Mi  =  I  and  -  0^ . 


By  introducing  the  transformation  q  =  ii)  Eq.  (4.2.1)  becomes 

q  +  O'^q  =  *^BpU 


(4.2.28) 


y  =  ^y*r\ 


(4.2.2b) 


A  state  space  model  of  the  system  in  Eq.  (4.2.2)  can  now  be 
constrncted.  Let 


q 

L’' 


Then  Eq,  (4.2.2)  can  be  written  in  state  space  form. 


q 

1 

O 

M 

_ 1 

q 

0 

— 

+ 

•  e 

q 

1 

K> 

O 

J _ _ 

q 

(4.2.3a) 


y  =  10  Cl] 


[3J  • 


(4.2.3b) 


If  N  is  positive  definite  then  another  state  space  description  can  be 
obtained  by  multiplying  through  in  Eq.  (4.2.1)  by  N~^ .  Letting 


yields 
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“  - 

4 

0  I 

q 

+ 

0 

q 

-m"1k  0 

q 

BT^Bp 

y 


[0 


(4.2.4b) 


The  latter  state  space  description  was  used  in  the  design  process 

T  *T  T 

because  the  original  coordinates  (q  .q  )  are  retained  in  the  state  space 
descr ipt ion. 


For  ease  in  discussion,  the  following  equations  are  introduced,  with 
the  obvious  equivalence  with  Eq.  (4,2.4) 


z 


Ax  +  Bu 


(4.2.5a) 


y  =  Cx  .  (4,2.5b) 

4.3 .  A  Physical  Problem  Description 

The  structure  analyzed  herein  has  been  proposed  by  Charles  Stark  Draper 
Labs  and  can  be  seen  in  Figure  4.1.  All  of  the  numerical  values  are 
summarized  in  Appendix  B.  This  same  structure  has  been  analyzed  in  [19] 
using  positivity  concepts. 

The  tetrahedral  apex  represents  the  antennae  feed,  with  members  1-6 
forming  the  support  structure  and  bi-peds  7-8,  9-10,  and  11-12  being 
snpports/controls  which  are  fixed  to  an  inertially  stabilized  (assumed) 
antenna  dish.  The  physical  nature  of  the  problem  allows  sensors  to  be 
placed  only  on  the  bi-peds;  no  sensors  can  be  placed  on  the  feed  itself. 


(?)  Leg  1  Station  1 


Member 

Number 


FP-  F867 


Figure  4.1:  Draper  Tetrahedral  Truss 
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The  actuators  control  the  elongation  and  contraction  of  the  bi-peds.  Only 
velocity  information  has  been  used  in  the  system  model;  displacement  of  the 
bi-peds  is  not  sensed  or  controlled. 

Because  uncertainty  in  the  model  parameters  exists,  a  perturbed  system 
must  be  studied.  A  perturbed  system  is  siaiply  different  M  and  K  matrices. 
For  the  design  to  be  acceptable,  the  objective  must  be  satisfied  for  both 
the  nominal  and  perturbed  systems. 

The  objective  is  to  damp  the  x  and  y  deflection  of  the  feed,  node  1,  to 
less  than  .0004  and  .00025  units,  respectively,  in  20  seconds.  Thus,  the 
problem  is  to  control  the  apex  in  the  presence  of  modelling  uncertainty  and 
without  directly  controlling  or  measuring  its  motion  or  position. 

The  Design  Procedure 

The  design  process  began  by  verifying  the  observability  and 
controllability  of  both  the  nominal  and  perturbed  systems  using  the  chained 
aggregation  algorithm.  The  eigenvalues  of  the  two  systems  were  also 
computed  and  can  be  found  in  Table  4.1. 

The  next  step  was  to  obtain  a  reduced  order  model  which  accurately 
represented  the  overall  system  in  general  and  included  the  desired  x  and  y 
deflection  modes  of  node  1  in  particular. 

The  original  output  structure,  the  C  matrix,  of  the  system  contained 
only  velocity  measurements  of  the  three  bi-peds  and  by  aggregating  the 
system  with  respect  to  this  information,  the  x  and  y  deflections  of  node  1 
would  be  forced  into  the  residual  subsystem,  an  undesirable  result.  To 
circumvent  this  problem  a  C  matrix  containing  both  the  controlled  outputs. 
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Table  4.1;  Nominal  and  Perturbed  System  Eigenvalues 


SYSTEM  EIGENVALUES 


NOMINAL 

0  +  jl2.92 

0  +  jlO.28 

0  +  j9.25 

0  +  j8.54 

0  +  j4 .76 

0  +  j4.66 

0  +  j4.20 

0  ±  j3.40 

0  +  j2.96 

0  +  j2.89 

0  ±  jl.66 

0  ±  jl.34 

PERTURBED 

0  +  jl3.97 

0  +  jlO.92 

0  +  jl0.30 

0  +  j8.94 

0  ±  j5.71 

0  +  j5.68 

0  +  j5.15 

0  +  j3.85 

0  +  j3.56 

0  +  j2.96 

0  +  jl.47 

0  +  jl.l7 

the  z  and  y  displacement  of  node  1,  and  the  measured  outputs,  the  original, 
nominal  C  matrix,  was  constructed. 

Using  this  C  matrix  together  with  the  nominal  system  matrix  A,  four 
steps  of  the  chained  aggregation  algorithm  were  implemented  with  the  system 
not  aggregating.  The  transformation  matrix  resulting  from  this  aggregation 
process  was  then  used  to  transform  the  nominal  system,  (A,B,C) ,  and  the 
perturbed  system,  (A^,bP,C^).  The  following  structure  for  both  systems  was 
ident if ied 


^12] 

, 

BQ 

L _ 

A«  „ 

z 

0 

22 

r  1 

•  _ 

u 


(4.4.1a) 
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y  =  [c, 


“'L') 


(4.4.1b) 


where  G  z^  £  R®,  u  €  ,  y  €  R^,  and  the  oatrices  are  partitioned 
accordingly. 

The  eigenvalues  of  the  and  ^22  sob“*lrices,  subsequently  referred 
to  as  the  aggregate  and  residual  subsystems,  respectively,  were  computed  and 
are  in  Table  4.2.  Two  innaediate  conclusions  are  obtained  by  comparing  the 
eigenvalues  in  Table  4.1  with  those  in  Table  4.2.  First,  the  two  sets  ! 
eigenvalues  are  both  purely  imaginary  and  nearly  equal.  Second,  the  lower 
frequencies  were  not  necessarily  placed  in  the  aggregate  subsystem,  but 


Table  4.2;  Nominal  and  Perturbed  Subsystem  Eigenvalues 
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instead  the  strongly  observable  modes  have  been  forced  into  the  aggregate. 

This  differs  from  a  modal  decomposition  of  the  system  where  the  lower 
frequency  modes  are  retained  to  form  a  reduced  order  model  and  the  remaining 
modes  are  discarded  [24].  With  the  chained  aggregation  procedure  the 
important  modes  are  those  modes  which  influence  the  observable  modes  of  the 
system  the  most.  It  is  not  necessarily  concerned  with  the  relative 
frequencies  of  the  modes. 

The  effect  of  the  neglected  modes,  both  residual  and  infinite 
dimensional,  on  the  reduced  order  model  has  been  referred  to  in  the 
literature  as  controller  and  observer  spillover  [17,18],  This  refers  to  how 
the  unmodelled  modes  are  affected  by  the  input  and  how  they  affect  the 
output,  respectively. 

An  initial  check  to  see  how  much  coupling  existed  from  the  residual 
subsystem  into  the  aggregate  subsystem  was  performed  by  comparing  the  sizes 
of  the  submatrices  involved.  The  size  used  for  each  submatrix  was  the 
matrix  2-norm,  i.e.,  the  largest  singular  value.  The  norms  computed  were 

I  I  I  1^  =  81  .92  ,  I  I  I  *2  " 

II  A21  II2  =  94.40  ,  II  A22  II2  =  148.37  . 

and 

II  AP^  I =  94,13  . 


II  A^2  **2  = 
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K  A?!  II2  =  97.88  .  II  a|2  ^^2  “  178.18  . 


The  differences  in  magnitude  suggested  the  coupling  was  weak  and  could  be 
neglected  [25] . 


Further  analysis  supporting  this  claim  was  obtained  from  a  geometric 
setting,  following  the  work  in  [61  on  near  unobservability.  Intuitively,  if 
the  residual  subsystem  was  nearly  unobservable  in  the  aggregate  subsystem, 
then  the  coupling  of  the  residual  into  the  aggregate  would  be  considered 
weak,  i.e.,  the  ^12  submatrix  would  have  little  influence  on  the  aggregate 
states.  As  shown  later  this  was  made  more  rigorous  by  showing  that  the 
subspace 


was  near  the  snbspace 

To  1 


This  suggested  that  the  residual  was  nearly  unobservable  by  the  aggregate, 
so  the  snbmatrix  was  neglected  in  the  analysis. 

To  strengthen  the  concept  of  one  snbspace  being  near  another,  canonical 
angles  are  defined. 


Definition  4.4..1  [6.261  :  Let  IT  and  V  be  subspac. .  jf  R°  with  orthonormal 

T 

bases  U  and  V,  respectively.  Let  or-  be  the  singular  values  of  U  V.  Then 
the  canonical  angles  between  U  and  T  are  the  numbers 


I 
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0^  =  Cos"^  □ 

Referring  to  Figure  2.1  a  basis  for  D  and  a  basis  for  V  was  obtained  by 

T  XX 

performing  singular  value  decomposition  on  both  matrices 
T  T 

[0  '^22^  '  using  the  generated  submatriz  for  each  basis. 

The  canonical  angles  between  snbspaces  D  and  V  were  then  obtained  using 
Definition  (4.4.1)  and  can  be  found  in  Table  4.3  for  both  the  nominal  and 
perturbed  systems.  In  general,  the  canonical  angles  were  small  and  the  two 
subspaces  were  considered  near  each  other. 

Another  motivation  for  neglecting  the  coupling  due  to  the  submatriz 

is  because  of  the  structure  of  the  transformed  C  matriz  and  the  fact  that 
output  feedback  is  being  used,  i.e.,  the  A^^  submatriz  will  not  be  affected 
by  the  feedback. 


Table  4.3;  Canonical  Angles  Between  D  and  ▼  For 
The  Nominal  and  Perturbed  Systems 


CANWJICAL  ANGLES 


0 

0 

«  «  ® 

0 

NOMINAL 

0.0 

0.0 

0.0 

7  .3 

20.5  * 

29.9  * 

34.6  " 

0 

42.4 

PERTURBED 

0.0  " 

• 

0.0 

12.6  * 

0 

13  .0 

0 

0 

0 

21.5 

29.9 

34.6 
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By  neglecting  the  coupling  end  noting  that  the  residual  subsysteia  was 
stable,  the  aggregate  subsystem  yielded  the  following  reduced  order  model 
which  was  used  in  the  design  process. 

*a  =  +  B^u  (4.4.2a) 

y  =  .  (4.4.2b) 

This  model  as  well  as  the  reduced  order  perturbed  model  were  verified  to  be 
controllable  and  observable  by  using  the  chained  aggregation  algorithm. 

Because  all  of  the  eigenvalues  in  the  aggregate  were  imaginary,  some 
type  of  damping  had  to  be  introduced.  The  structure  of  the  transformed 
nominal  input  and  nominal  output  matrices  suggested  that  output  feedback, 
u  -  -  Ky,  be  used  to  introduce  the  desired  damping.  The  gain  matrix,  K,  was 
obtained  using  the  procedure  outlined  in  [27,28];  more  background  material 
is  given  in  [29,30] .  This  design  approach  begins  by  solving  a  reference 
optimal  state-feedback  linear  quadratic  regulator  problem.  The  eigenvectors 
associated  with  the  closed-loop  system  matrix  are  then  determined.  If  there 
are  r  system  outputs,  then  r  of  these  eigenvectors  are  retained  in  the 
reduced-order  output  feedback  problem,  i.e.,  an  r-dimensional  eigenspace  of 
the  reference  problem  is  retained. 

To  achieve  the  desired  objective  for  the  x  and  y  deflections  of  node 
one,  the  appropriate  Q  and  R  matrices  had  to  be  selected.  The  following 
analysis  led  to  the  desired  weighting  matrices.  The  input  weightings  were 
selected  to  be  equal  and  large  so  that  the  input  energy  would  not  be 


excessive.  The  weighting  on  the  states  was  not  as  straightforward. 


After 
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the  necessary  transformations  were  carried  out  to  place  the  system  in  the 
correct  structure,  the  x  and  y  deflections  were  composed  largely  of  states 
15  and  16  with  states  13  and  14  corresponding  to  their  respective 
velocities.  Because  the  z  and  y  deflections  as  well  as  their  velocities 
were  nearly  uncontrollable  in  the  reduced-order  model,  less  weighting  than 
might  be  expected  on  these  states  had  to  be  used.  Looking  at  the  feedback 
structure,  states  1-6  were  weighted  very  heavily  and  states  7-11  were  hardly 
weighted,  by  comparison.  State  12  did  require  more  damping,  so  it  was 
weighted  appropriately.  These  conclusions  were  drawn  after  preliminary 
trials  and  observing  how  the  various  states  related  to  one  another  and  how 
their  magnitudes  compared  with  the  other  states.  The  resulting  x  and  y 
deflections  for  both  reduced  order  models  can  be  seen  in  Figures  4.2  and 
4.3. 

To  verify  that  the  above  gain  matrix,  K,  resulted  in  a  desirable  design 
for  the  full  order  systems,  it  was  applied  to  both  the  nominal  and  perturbed 
full  order  models.  The  resulting  x  and  y  deflections  can  be  seen  in  Figures 
4.4  and  4.5. 

By  examining  the  trajectories  of  the  full  order  systems  one  can  see  the 
slight  error  introduced  by  using  the  reduced-order  model  during  the  design 
procedure.  The  x  deflections  are  seen  to  meet  the  required  objective,  but 
the  y  deflections  require  28  and  34  seconds  to  meet  the  stricter  objective 
in  the  nominal  and  perturbed  systems,  respect ivelv . 

The  use  of  suboptimal  dynamic  output  feedback  [27,28,29]  could  possibly 
increase  the  flexibility  and  robustness  of  the  solution  and  should  be 
investigated  in  future  work. 


X10 


igure  4.3  Trajectories  for  the  perturbed  reduced-order  model 


CHAPTER  5 


CONCLDSION 


A  computer  implementation  of  chained  aggregation  and  modified  chained 
aggregation  using  orthogonal  transformation  matrices  has  been  presented.  To 
obtain  the  orthogonal  matrices  the  singular  value  decomposition  has  been 
incorporated  into  the  algorithm  whenever  a  state  space  transformation  must 
be  performed.  Because  orthogonal  matrices  have  been  used,  the  problem  sen¬ 
sitivity  will  usually  not  be  altered. 

In  Chapter  4  the  algorithm  has  been  used  in  the  design  of  a  control  for 
a  particular  large  space  structure  prototype.  The  reduced-order  model  sug¬ 
gested  by  using  the  chained  aggregation  algorithm  contains  the  desired 
information  structure,  the  x  and  y  deflections  of  node  one,  and  can  be  used 
during  the  design  process  and,  or  in  the  physical  implementation  of  the  con¬ 
trol,  if  a  dynamical  feedback  is  to  be  used. 

Future  research  may  now  focus  on  the  numerical  solution  of  practical 
problems  using  the  design  schemes  presented  in  past  Generalized  Hessenberg 
Representation  articles. 


APPENDIX  A 


SUPPORTING  SOFTWARE 


The  two  algorithms  developed  in  Chapters  2  and  3  have  been  incorporated 
into  a  single  computer  program.  The  program  has  been  written  in  FORTRAN 
with  the  singular  value  decomposition  subroutine  taken  from  LINPACK  [31].  A 
flowchart  of  the  software  follows. 


Flowchart  continued  on  next  page. 


Flowchart  continued  on  next  page 
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System  will  not  aggregate. 


No  C  matrix. 
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3 


II 

[j 

ij  An  explanation  of  the  variable  names  follows; 


a  =  The  number  of  linearly  ndependent  columns  in  the  C  and 
A^2  submatrices.  It  is  redefined  after  every  column 
compression  which  is  performed  on  these  submatrices. 

P  =  The  number  of  linearly  independent  rows  in  the  B. 

submatrix.  It  is  redefined  after  every  compression  of 
the  rows  of  Bj . 

T  =  The  sum  of  the  a's  as  they  are  generated. 

p  =  The  size  of  the  current  residual  subsystem. 

n  =  The  number  of  zero  rows  identified  in  the  top  of  the  B 

matrix . 


The  equivalence  between  the  greek  symbols  used  herein  and  the  variable 
names  used  in  the  actual  computer  program  are  as  follows. 

T  =  NIDE 
^  =  NOZERO 

a  =  NZ 
P  >=  NZB 
p  =  SIZA 

During  the  design  process  in  Chapter  4  additional  software  was 
developed.  Many  of  the  routine  linear  matrix  manipulations  were  performed 
using  LAS  132].  The  software  required  to  compute  the  gain  matrix,  K,  was  a 
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mixture  of  LAS  code  and  FORTRAN  software.  The  IMSL  subroutine  EIGRF  was 
used  to  obtain  the  eigenvalues  and  eigenvectors  of  the  closed-loop  matrix. 
A  flowchart  description  of  the  steps  involved  follows; 
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Select  r  of  these  eigenvectors,  form  W. 
(for  complex  conjngate  pairs 
the  real  and  imaginary 
parts  make  np  2  of  the 
columns) 


Flowchart  continued  on  next  page. 


Acceptable 


GO  TO  START 


GO  TO  A 


Software  independent  of  LAS  was  written  to  generate  the  system 

trajectories  using  the  initial  condition,  the  state  transition  matrix,  the 

total  time,  and  the  step  size.  The  state  transition  matrix  was  obtained 
At 

using  the  e  operator  in  LAS.  A  graphing  routine  which  allows  the  plotting 
of  any  of  the  resulting  system  trajectories,  as  well  as  the  x  and  y 
deflections  of  node  one  was  also  developed. 
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APPENDIX  B 

NUMERICAL  VALDES  ASSOCIATED  WITH  THE  LARGE  SPACE  STRUCTURE 


The  niunerical  values  associated  with  Eq.  (4.2.1)  for  both  the  nominal 
and  perturbed  systems  are  given  in  Tables  B.l ,B.2,B.3 ,B.4  and  Tables 
B.S,B.6i  respectively.  Note  that  matrices  Bp  and  C^  are  identical  for  both 
systems.  If  the  state  space  representation  in  Eq.  (4.2.3)  is  desired,  the 
required  d  matrices  are  given  in  Tables  B.7  and  B.8.  The  relationship 

between  the  coordinates  (q  »  q  )^  ^nd  the  nodes  of  Figure  4.1  can  be  found 
in  Table  B.9.  The  transformation  matrix  used  to  generate  Eq.  (4.4.1)  for 
both  the  nominal  and  perturbed  systems  is  given  in  Table  B.IO.  All  of  the 
reduced-order  system  matrices  (Aj^.Bj.Cj)  and  (A^j.B^.CP)  are  contained  in 
Tables  B. 11 , B. 12 , B. 13  and  Tables  B.14,B.1S,  and  B,16,  respectively.  To 
obtain  the  system  used  during  the  design  of  the  K  feedback  matrix,  the 
transformation  matrix  in  Table  B.17  was  used.  Tables  B.18  and  B.19  contain 
the  final  Q  and  R  matrices  used  to  generate  the  K  feedback  matrix  which  is 
given  in  Table  B.20. 
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Table  B.l:  Mass  Matrix  for  the  Nominal  System 


2  0  0  0  0  0 
0  2  0  0  0  0 
0  0  2  0  0  0 
0  0  0  2  0  0 
0  0  0  0  2  0 
0  0  0  0  0  2 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 


0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0  0 
0  0  0  0  0 
0  0  0  0  0 
2  0  0  0  0 
0  2  0  0  0 
0  0  2  0  0 
0  0  0  2  0 
0  0  0  0  2 
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Table  B.3:  The  Bp  Matrix 


0000 

.0000 

.0000 

.0000 

.0000 

.0000 

0000 

.0000 

.0000 

.0000 

.0000 

.0000 

0000 

.0000 

.0000 

.0000 

.0000 

.0000 

3535 

.3536 

.0000 

.0000 

.0000 

.0000 

6124 

-.6123 

.0000 

.0000 

.0000 

.0000 

7071 

-.7071 

.0000 

.0000 

.0000 

.0000 

0000 

.0000 

-.3536 

.3535 

.0000 

.0000 

0000 

.0000 

-.6123 

.6124 

.0000 

.0000 

0000 

.0000 

-.7071 

-.7071 

.0000 

.0000 

0000 

.0000 

.0000 

.0000 

.7071 

-.7071 

0000 

.0000 

.0000 

.0000 

.0000 

.0000 

0000 

.0000 

.0000 

.0000 

-.7071 

-.7071 

Table  B.9:  Relationship  Between  the  Coordinates  and  the  Nodes 


<11 

( qj)  = 

Displacement 

(Velocity) 

of 

node 

1 

in 

the 

negative 

X 

direction. 

‘^2 

(qj)  = 

Displacement 

(Velocity) 

of 

node 

1 

in 

the 

negative 

y 

direct  ion. 

^3 

(qs)  = 

Displacement 

(Velocity) 

of 

node 

1 

in 

the 

negative 

d irect ion. 

(q4)  = 

Displacement 

(Velocity) 

of 

node 

2 

in 

the 

nega  t ive 

X 

direction. 

“If 

(\s)  5 

Displacement 

( Veloc ity) 

of 

node 

2 

in 

the 

nega  t ive 

y 

direct  ion. 

(q<;)  = 

Displacement 

(Velocity) 

of 

node 

2 

in 

the 

negative 

z 

direction. 

^7)  S 

Displacement 

(Velocity) 

of 

node 

3 

in 

the 

negative 

X 

direct  ion. 

'^8 

( qg)  = 

Displacement 

(Velocity) 

of 

node 

3 

in 

the 

negative 

y 

direct  ion. 

<  q9)  S 

Displacement 

(Velocity) 

of 

node 

3 

in 

the 

negative 

z 

direct  ion. 

1  **110^ 

=  Displacement  (Velocity)  ' 

of  node 

4 

in  the  negative 

z  direction. 

‘^11 

.  ^^11^ 

—  Displacement  (Velocity)  1 

of  node 

4 

in  the  negative 

y  direction. 

‘^12 

;  ^<112^ 

—  Displacement  (Velocity)  ' 

of  node 

4 

in  the  negative 

z  direction. 
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Table  B.IO:  Transformat 

T 

ion  matrix  which 
=  t  Tj.  T3  1 

generates  Eq. 

(4.4.1) 

r 

1 .0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

1.0000 

.0000 

.0000 

.0000 

.0000 

,0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

,0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

,0000 

.0000 

.0000 

,0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

,0000 

.0000 

.0000 

.0000 

.0000 

,0000 

,0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.3104 

-.4680 

-.2701 

-.4736 

.0660 

.4352 

.0000 

.0000 

-.0002 

-.3374 

.5846 

.0003 

.5133 

.0934 

.0000 

.0000 

.4170 

,0379 

.0222 

-.6324 

-.0673 

-.4554 

.0000 

.0000 

.1492 

.0171 

.0101 

.1074 

.0730 

.0017 

.0000 

.0000 

-.0001 

-.0236 

.0406 

.0000 

-.0325 

.0460 

.0000 

.0000 

.0000 

-.2887 

.5000 

,0000 

-.2887 

-.5000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

,0000 

.0000 

.0000 

.0000 

.0000 

.5646 

-.4162 

-.2401 

.5909 

-.1369 

-.1567 

.0000 

.0000 

-.0003 

-.2285 

.3958 

0002 

-.3900 

.4890 

.0000 

.0000 

,0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

-.0044  - 

.0058 

-.0034 

.0002 

.0076 

-.0036 

.0000 

.0000 

-.2904  - 

.3933 

-.2272 

.0200 

.5403 

-.2494 

.0000 

.0000 

-.5517  - 

.4487 

-.2591  - 

.1211  - 

-.4185 

-.1283 
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Table  B.IO:  (Contiuned) 


.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.000^' 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

-.0030 

.4097 

-.1606 

-.0030 

-.0049 

-.3374 

-.3979 

.0049 

-.0118 

-.4283 

.1691 

-.0118 

-.6905 

.0381 

.0624 

-.6916 

.7046 

-.0237 

.0513 

-.7036 

.0003 

.5773 

-.0001 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0851 

-.2043 

-.0398 

.0852 

-.0653 

-.2285 

.5824 

.0652 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0002 

.0007 

.0084 

.0002 

.0200 

.0543 

.5924 

.0203 

-•1210  -.3203  -.2984  -.1210 


.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

-.2192 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

,0000 

.0000 

.0000 

-.4494 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

,0000 

.0000 

.0000  ; 

.0000 

.0000 

.0000 

( 

,0000  { 

0000 

.0000 

.0000 

1 

.0000  j 

-.9996  -.0288  .0000  .0000 

.0288  -.9996  .0000  .0000 

.0000  .0000  .0000  .0000 

.0000  .0000.  .0000  .0000 


.0000 

.0000 

-.6203 

.6515 

.0000 

.0000 

.0002 

.0001 

.0000 

.0000 

-.6875 

-.3336 

.0000 

.0000 

.3776 

.4628 

.0000 

.0000 

-.0055 

-.0067 

,0000 

.0000 

-.0001 

-.0002 

J 
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Table  B.IO:  (Continued) 


p  « 

.0000 

.0000 

.0000 

.0000  .0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000  .0000 

.0000 

.0000 

.0000 

.0000 

.0000 

-.2188 

-.3792  .8991 

.0000 

.0000 

.0000 

.3197 

-.8988 

.0000 

.0000  .0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000  .0000 

.7071 

.0000 

-.1011 

.0000 

.0000 

.0000 

.0000  .0000 

-.7071 

.0000 

-.1071 

.1783 

.4384 

.0000 

.0000  .0000 

.0000 

.0000 

.0000 

.0000 

.0000 

-.4496 

-.7786  -.4378 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000  .0000 

.0000 

,0000 

.0000 

1  .0000 

.0000 

.0000 

.0000  .0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000  .0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000  .0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

•0000  .0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000  .0000 

.0000 

,0000 

.0000 

.0000 

.0000 

.0000 

.0000  .0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000  .0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000  .0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000  .0000 

.0000 

.0000 

.0000 

.3761 

.0001 

.1360  - 

•0785  .0000 

.0000 

.1569 

.0000 

.0000 

.0000 

.6124  - 

.3536  -.0001 

.0000  - 

.7070 

.0000 

-.1926 

.0000 

.3769 

.2176  .0000 

.0000  -, 

.43  55 

1 

,0000  1 

.2672 

0001 

4630 

2673  .0001 

0000  -. 

5346 

1 

J 

0000  1 

-.0039 

0000 

0067 

0039  .0000 

0000 

0077 

1 

1 

0000  ! 

-.0001  .0000 

0002 

0001  .0000 

0000 

0002 

1 

0000  i 
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Table  B.ll:  Nominal  rednced-order  system  matrix 
A  —  I  Aj ,  A2  ] 


.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

-1.83 

1.06 

.00 

.00 

.00 

.00 

.00 

.00 

-18.34 

-10.59 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

1.67 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

1.67 

.00 

.00 

.00 

.00 

.00 

.00 

8.95 

5.17 

.00 

.00 

.00 

.00 

.00 

.00 

.89 

-.52 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

-.18 

.43 

.29 

-.29 

-.21 

-.37 

.00 

.00 

-.15 

.30 

-.24 

.24 

.61 

-.30 

.00 

.00 

.18 

.57 

-.29 

.31 

-.28 

.39 

.00 

.00 

-.64 

-.10 

.52 

.46 

.05 

.27 

.00 

.00 

.62 

.02 

.48 

.51 

.04 

-.33 

.00 

.00 

.25 

.25 

.41 

-.41  • 

.52 

.52 

13.93 

6.67 

.00 

.00 

.00 

.00 

.00 

.00 

6.10 

6.06 

.00 

.00 

.00 

.00 

.00 

.00 
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Table  B.ll:  (Continued) 


.00 

.00 

.00 

.00 

.00 

.00 

-1.00 

.03 

.00 

.00 

.00 

.00 

.00 

.00 

-.03 

-1.00 

17.19 

12.86 

-10.53 

13.76 

-13.33 

-4.47 

.00 

.00 

-42.45 

-23.86 

-40.56 

.87 

-.17 

-4.48 

.00 

.00 

-11 .11 

19.99 

18.54 

-10.96 

-10.40 

-7.22 

.00 

.00 

26.20 

-18.51 

-21.41 

-10.50 

-10.67 

7.22 

.00 

.00 

10.10 

-48.97 

19.77 

-.43 

-.34 

-9.17 

.00 

.00 

33.66 

24.04 

-26.28 

-6.35 

6.71 

-9.17 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

23.77 

-.57 

19.98 

-1.26 

1.44 

.00 

.00 

.00 

12.83 

1.02 

10.78 

-.68 

-2.66 

.00 

.00 

.00 
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Table  B.12: 

Input  matrix  for 

the  nominal  reduced 

-order  model 

r 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

-.1631 

-.4726 

.0000 

-0000 

.4728 

.1628 

.0000 

.0000 

.0000 

,0000 

.0000 

.0000 

.0000 

.0000 

.5000 

.0000 

.0000 

.0000 

.0000 

.0000 

,0000 

.5000 

.1627 

-.4728 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.4726 

-.1631 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

1 

.0000  1 

.0000 

.0000 

.0000 

,0000 

.0000 

.0000  j 

.0000 

.0000 

.0000 

.0000 

.0000 

1 

. 0000  1 

.0000 

.0000 

,0000 

.0000 

.0000 

t 

.0000  1 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000  j 
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Table  B.13: 


Output  matrix  for 


the  nominal  reduced-order 


.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

-.3262 

.0000 

.0000 

-.9453 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

0000 

.0000 

.0000 

0000 

.0000 

.0000 

.0000  .0000  .0000 
.0000  .0000  .0000 

.0000  .0000  ,0000 


.9455 

.0000 

.0000 

.3255 

.3255 

.0000 

.0000 

-.9455 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000  1.0000 

.0000 

.0000 

.0000 

.0000 

1.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

model 


.0000 

.0000 

.9453 

-.3262 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

0000 

0000 


Table  B.14; 

Perturbed 

ap  = 

redaced- 
[  A2 

-order 

] 

system 

matrix 

.00 

.00 

.00 

.00 

.00 

.00 

,00 

.00 

.00 

.00 

.00 

.00 

.00 

,00 

.00 

.00 

-2.75 

1.59 

.00 

.00 

.00 

.00 

.00 

.00 

-22.01 

-12.71 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

2.50 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

2.50 

.00 

.00 

.00 

.00 

.00 

.00 

10.74 

6.20 

.00 

.00 

.00 

.00 

.00 

.00 

1.34 

-.77 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

-.18 

.43 

.29 

-.29 

-.21 

-.37 

.00 

.00 

-.15 

.30 

.24 

.24 

.61 

-.30 

.00 

.00 

.18 

.57 

.29 

.31 

-.28 

.39 

.00 

.00 

-.64 

-.10 

.52 

.46 

,05 

.27 

.00 

.00 

.62 

.02 

.48 

.51 

.04 

-.33 

.00 

.00 

.25 

.25 

.41 

-.41 

.52 

.52 

8.54 

3.91 

.00 

.00 

.00 

.00 

.00 

.00 

3.55 

3.95 

.00 

.00 

.00 

.00 

.00 

.00 

'  1! 

J 
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Table  B.14:  (Continued) 


.00 

.00 

.00 

.00 

.00 

.00 

-1.00 

.03 

.00 

.00 

.00 

.00 

.00 

.00 

-.03 

-1.00 

21.77 

16.49 

-13.08 

20.69 

-19.97 

-6.71 

.00 

.00 

-53.20 

-30.20 

-51.69 

1.56 

-.31 

-6.71 

.00 

.00 

-34.35 

25.48 

23.42 

-16.51 

-15.55 

-10.83 

.00 

.00 

32.83 

-23 .27 

-27.74 

-15.68 

-16.06 

10.83 

.00 

.00 

25.94 

-61.98 

25.19 

-.76 

-.63 

-13.76 

.00 

.00 

42.25 

30.32 

-33.84 

-9.43 

10.12 

-13.76 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

.00 

14.20 

-.42 

11.82 

-1.03 

1.08 

.00 

.00 

.00 

7.66 

.77 

6.38 

-.56 

-2.00 

.00 

.00 

.00 

72 


B.15:  Input  matrix  for  the  perturbed  reduced-order  model 


.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.4728 

.1628 

.0000 

.0000 

.0000 

.0000 

.1627 

-.4728 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.1631 

-.4726 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.4726 

-.1631 

.0000 

.0000 

.0000 

.0000 

,0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

,0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

,0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.5000 

.0000 

.0000 

.5000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 
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Table  B.16:  Output  matrix  for  the  perturbed  reduced-order  model 


CP 

=  f  s 

.  C2  ] 

0000 

.0000 

.0000 

.9455 

.0000 

.0000 

.3255 

.0000 

0000 

.0000 

.0000 

.3255 

.0000 

.0000 

-.9455 

.0000 

0000 

.0000 

-.3262 

.0000 

.0000 

.0000 

.0000 

.9453 

0000 

.0000 

-.9453 

.0000 

.0000 

.0000 

.0000 

-.3262 

0000 

.0000 

.0000 

.0000 

1.0000 

.0000 

.0000 

.0000 

0000 

.0000 

.0000 

.0000 

.0000 

1.0000 

.0000 

.0000 

0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 
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Table  B.17.  Transformation  nsea  duri 

T  =  [  ■ 


r 


.0000 

.0000 

.0000 

.9455 

.0000 

.0000 

.0000 

.3255 

.0000 

.0000 

-.3262 

.0000 

.0000 

.0000 

-.9453 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.2198 

.1269 

.0000 

.0000 

- . 0006 

.0008 

.0000 

.0000 

-.1990 

-.1149 

.0000 

.0000 

.0950 

.0549 

.0000 

.0000 

-.0678 

.1173 

.0000 

.0000 

-.0001 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

-.9082  ■ 

-.2308 

.0000 

.0000 

-.2715 

.9490 

.0000 

.0000 

the  des 
.  T2  J 

ign  of 

the  feedback  ma 

.0000 

.0000 

.3255 

.0000 

.0000 

.0000 

-.9455 

.0000 

.0000 

.0000 

.0000 

.9453 

.0000 

.0000 

.0000 

-.3262 

1.0000 

.0000 

.0000 

.0000 

.0000 

1.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

Table  B.17:  (Continued) 


.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.9549 

-.0001 

.1541 

-.0030 

.0002 

1.0000 

-.0001 

.0000 

.2143 

-.0002 

-.9493 

.0039 

.0264 

.0000 

-.0366 

-.9929 

.0000 

-.0001 

.0000 

-.0001 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.1970 

-.0003 

.2623 

-.1146 

.0523 

-.0009 

-.0696 

.0304 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

,0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0001 

.0000 

.0000 

.0000 

.9908 

.0000 

.0000 

.0000 

.0000 

-1.0000 

.0000 

.0000 

.0000 

,0000 

-.8801 

-.4748 

.0000 

.0000 

-.4748 

.8801 

.0348 

.0001 

.0000 

.0000 

.1310 

.0001 

.0000 

.0000 
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Table  B.18:  State  weighting  matrix  U 

10^  000000000000000 
0  10^  00000000000000 
00  10^  0000000000000 
000  10^  000000000000 
0000  10*^  00000000000 
0  0  0  0  0  10^  0000000000 
000000  10^  000000000 
0000000  10^  00000000 
00  )  00000  10^  0000000 
000000000  10^  000000 
0000000000  10^  000  0 
00000000000  10'*  0  0  0  0 
000000000000  lO'*  000 
0000000000000  10“*  00 
00000000000000  10'*  0 
000000000000000  10'* 
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Table  B.19:  Input  weighting  aatria  R 


R  = 


10^  0 
0  10^ 
0  0 
0  0 
0  0 
0  0 


0  0 
0  0 
10^  0 
0  10^ 
0  0 
0  0 


0  0 
0  0 
0  0 
0  0 
10^  0 
0  10^ 


Table 

B.20: 

Output 

feedback  gain  matrix 

K 

31.6249 

.0024 

.0087 

-.0134 

-.0140 

.0086 

.0024  31 

.6249 

.0086 

-.0140 

-.0134 

.0087 

-.0039 

.0025 

31.6224 

.0006 

.0011 

-.0005 

.0125 

.0134 

-.0002 

31.6228 

.0006 

-.0013 

.0134 

.0125 

-.0013 

.0006 

31.6228 

-.0002 

-.0025  -.0039  -.0005 


0011 


0006  31.6224 
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APPENDIX  C 


PROGRAM  LISTING 


CHARACTER*! 

ANS  .  CFRMT 

AFRMT 

,  BFRMT 

CHARACTER*20 

FANAME  .  FCNAME 

FBNAME 

.  NAME 

CHARACrER*20 

INAME 

INTEGER 

I  .  J 

NCC 

,  NRC 

INTEGER 

NM  .  NMIN 

lERR 

.  NOZERO 

INTEGER 

FA12  ,  NIDE 

NZ 

.  SOS 

INTEGER 

NRB  .  NCB 

SIZA 

,  NZB 

INTEGER 

R  ,  ID 

JD 

,  NRGDMB 

INTEGER 

NMINB  .  FA12S 

NRA12S 

.  STEP 

INTEGER 

REAL  WORK(IOO) 

FCHND  ,  NCBTMP 

REAL  A(IOO.IOO)  . 

B( 100, 100) 

GdOO, 

100) 

REAL  C( 100, 100)  . 

A12(100.100) 

TA12d00.100) 

REAL  UdOO.lOO)  . 

VdOO.lOO) 

TdOO, 

100) 

REAL  E(IOO) 

GDMBdOO.lOO) 

TDMBdOO.lOO) 

REAL  A12B(100,100) 

.  A12Sd00d00) 

DMBdOO.lOO) 

REAL  NORM 

REAL  SIGMA(IOO) 

SIZCHX 

LOGICAL 

INFO 

VERSION  aggregate.!  CREATED  DECEMBER  22,1982  BY:  HAL  THARP  • 


TO  COMPILE  THIS  PROGRAM; 
ill  -u  -o  agg  aggregate.! 


VARIABLE  TABLE:  MAIN  PROGRAM  • 


A  :  CONTAINS  THE  A  MATRIX  IN  THE  STATE  SPACE  DESCRIPTION. 

ANS  :  A  CHARACTER  VARIABLE  WHICH  IS  USED  TO  CONTAIN  THE 

RESPONSE  OF  THE  USER  TO  A  QUESTION. 


A12  CONTAINS  THE  CURRENT  A12  SUBMATRIX. 

A12B  CONTAINS  THE  SUBMATRIX  RESULTING  FROM  TRANSFORMING 

A12  WITH  PART  OF  THE  U  MATRIX  GENERATED  DURING  THE 
ROW  COMPRESSION  OF  A  B  SUBMATRIX.  THIS  MATRIX  IS 
CONTAINED  IN  THE  RANGE  OF  THE  INPUT. 


A12S 


CONTAINS  THE  SUBMATRIX  RESULTING  FROM  TRANSFORMING 
A12  WITH  PART  OF  THE  U  MATRIX  GENERATED  DURING  THE 


ROir  COMPRESSION  OF  A  B  SOBMAIRIX.  THIS  MATRIX  IS 
NOT  CONTAINED  IN  THE  RANGE  OF  THE  INPUT. 


B 

C 

DMB 

E 

FA12 

FA12S 

FCHND 

F_NAME 

_FRMT 

G 

GDNB 

lERR 

INAME 

INFO 

I.J.I 

ID.JD 


CONTAINS  THE  B  MATRIX  IN  THE  STATE  SPACE  DESCRIPTION. 

CONTAINS  THE  C  MATRIX  IN  THE  STATE  SPACE  DESCRIPTION. 

A  DUPLICATE  OF  GDMB  WHICH  IS  PASSED  TO  THE  SSVDC 
SUBROUTINE.  THIS  MUST  BE  DONE  BECAUSE  THE  MATRIX 
RETURNED  FROM  SSVDC  IS  ALTERED  FROM  THE  ORIGINAL 
MATRIX  PASSED. 

VECTOR  THAT  ORDINARILY  CONTAINS  ZEROES.  IT  IS  USED 
IN  THE  SSVDC  ROUTINE.  FOR  MORE  INFORMATION  CONSULT 
LINPACK  USERS'  GUIDE.  DONGARRA.  et.al.  CHPT  11. 

A  FLAG  WHICH  IS  ZERO  IF  THE  A12  SUBMATRIX  IS  ZERO 
AND  ONE  OTHERWISE. 

A  FLAG  WHICH  IS  ZERO  IF  THE  A12S  SUBMATRIX  IS  ZERO 
AND  ONE  OTHERWISE. 

A  FLAG  WHICH  IS  SET  TO  1  WHEN  CHAINED  AGGREGATION  IS 
PERFORMED  AND  0  WHEN  MODIFIED  CHAINED  AGGREGATION  IS 
PERFORMED. 

CONTAINS  THE  NAME  OF  THE  FILE  WHERE  THE  (A.B.C)  MATRICES 
ARE  STORED. 

CONTAINS  THE  FORMAT  METHOD  USED  IN  THE  STORAGE  OF  THE 
PARTICULAR  (A.B.C)  MATRIX,  a. A  =  SPARSE  OR  b.B  = 
NORMAL. 

THE  ARRAY  USED  TO  STORE  THE  SUBMATRIX  GENERATED  DURING 
THE  TRANSFORMATION  OF  THE  B  MATRIX.  THIS  SUBMATRIX  IS 
THEN  PLACED  IN  THE  CORRECT  LOCATION  WITHIN  B. 

A  TEMPORARY  ARRAY  WHICH  CONTAINS  THE  SUBMATRIX  OF  THE 
B  MATRIX  WHOSE  ROWS  ARE  TO  BE  COMPRESSED. 

VARIABLE  CONTAINING  ZERO  WHEN  THE  SINGULAR  VALUES  HAVE 
BEEN  COMPUTED  CORRECTLY.  SEE  LINPACK  USERS'  GUIDE. 

CHPT.  11  WHEN  lERR  IS  NOT  EQUAL  TO  ZERO. 

CHjIRACTER  VARIABLE  WHICH  IS  THE  NAME  OF  THE  FILE 
CONTAINING  THE  STAGE  INFORMATION.  IT  IS  NOT  USED  IF 
INFO  IS  FALSE. 

A  LOGICAL  VARIABLE  WHICH  IS  TRUE  WHEN  THE  USER  REQUESTS 
THE  STAGE  INFORMATION  TO  BE  OUTPUT  INTO  A  FILE. 

INTEGER  VARIABLES  USED  WITHIN  DO  LOOPS  AS  THE  COUNTERS. 
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NAME 

NCBTMP 

NC_ 

NIDE 

NM 

NMIN 

NMINB 

NORM 

NOZERO 

NR_ 

NZ 

NZB 

SIGMA 

SIZA 

SIZCHK 


A  'TIARACTEE  VARIABLE  CSED  TO  IDENTIFY  THE  PARTICULAR 

Matrix  being  input  into  ibe  program, 

AN  integer  variable  USED  TO  STORE  THE  NUMBER  OF 
COLUMNS  IN  TBE  ORIGINAL  B  MAIRII  WHEN  CHAINED  AGGREGATION 
IS  PERFORMED. 

:  CONSTANT  CONTAINING  THE  COLUMN  DIMENSION  OF  THE  PARTICULAR 
MATRIX  (B.C.ETC.). 

KEEPS  TRACK  OF  THE  NUMBER  OF  NON-ZERO  SINGULAR  VALUES 
GFTffiRATED  BY  COMPRESSING  THE  C  MATRIX  AND  THE  A12S 
SUBMATRICES.  THIS  IS  THE  NUMBER  OF  IDENTITY  ELEMENTS 
NEEDED  IN  THE  UPPER  LEFT  HAND  CORNER  WHEN  UPDATING  THE 
TRANSFORMATION  MATRIX, 

CONSTANT  CONTAINING  THE  MAXIMUM  POSSIBLE  DIMENSION  OF  THE 
MATRICES  PASSED  TO  THE  SINGULAR  VALUE  DECOMPOSITION 
ROUTINE  (SSVDC),  CURRENTLY  SET  EQUAL  TO  100. 

AN  INTEGER  VARIABLE  CONTAINING  THE  MINIMUM  OF  THE  NUMBER 
OF  COLUMNS  OH  THE  NUMBER  OF  ROWS  IN  THE  C  MATRIX  OR  AN 
A12S  SUBMATRIX. 

AN  INTEGER  VARIABLE  CONTAINING  THE  MINIMUM  OF  THE  NUMBER 
OF  COLUMNS  OR  THE  NUMBER  OF  ROWS  OF  THE  B  SUBMATRIX  WHICH 
IS  PASSED  TO  SSVDC.  THIS  IS  THE  MAXIMUM  NUMBER  OF  NON¬ 
ZERO  SINGULAR  VALUES  POSSIBLE  FOR  THE  PARTICULAR  SUBMATRIX 

CONSTANT  CONTAINING  THE  NORM  OF  THE  ORIGINAL  A  MATRIX. 

KEEPS  TRACK  OF  THE  NUMBER  OF  ZERO  ROWS  IN  THE  INPUT 
MATRIX. 

CONSTANT  CONTAINING  THE  ROW  DIMENSION  OF  THE  PARTICULAR 
MATRIX  (B.C.ETC.). 

AN  INTEGER  CONTAINING  THE  NUMBER  OF  NON-ZERO  SINGULAR 
VALUES  RESULTING  FROM  A  DECOMPOSITION  OF  THE  C  MATRIX 
OR  AN  A12  SUBMATRIX. 

NUMBER  OF  NON-ZERO  SINGULAR  VALUES  (STAINED  AFTER  PASSING 
PORTIONS  OF  THE  B  MATRIX  TO  THE  SSVDC  ROUTINE. 

;  VECTOR  CONTAINING  THE  SINGULAR  VALDES  RESULTING  FROM  THE 
SSVDC  ROUTINE, 

CONTAINS  THE  DIMENSION  OF  THE  A  MATRIX.  INITIALLY  IT  IS 
SET  EQUAL  TO  SOS,  BUT  THEN  IS  UPDATED  TO  KEEP  TRACK  OF  THE 
DIMENSION  OF  THE  A22  SUBMATRIX. 

CONSTANT  CONTAINING  THE  PRODUCT  OF  NORM  AND  THE  MACHINE 
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EPSILON.  ELEMENTS  SMALLER  TBAN  IBIS  IN  THE  SSVDC  ROUTINE 
ARE  TREATED  AS  BEING  ZERO. 

SOS  CONSTANT  CONTAINING  THE  DIMENSION  OF  THE  STATE  (SIZE  OF 

STATE) . 

STEP  KEEPS  TRACK  OF  THE  NUMBER  OF  MODIFIED  OR  CHAINED  STEPS 

PERFORMED. 

T  ARRAY  CONTAINING  THE  TRANSFORMATION  MATRIX  WHICH  IS 

CONTINUALLY  UPDATED  DURING  THE  ALGORITHM. 

TA12  TEMPORARY  STORAGE  ARRAY  WHICH  IS  USED  DURING  THE 

COMPUTATION  OF  A12.  SINCE  A12  IS  COMPUTED  IN  TWO  STAGES. 

TDMB  TEMPORARY  STORAGE  ARRAY  WHICH  IS  USED  DURING  THE  UPDATING 

OF  THE  TRANFORMATION  MATRTIX,  BECAUSE  THE  UPDATE  MUST  BE 
DONE  IN  TWO  STAGES. 

U  ;  ORTHOGONAL  MATRIX  WHOSE  COLUMNS  CONTAIN  THE  LEFT  SINGULAR 

VECTORS  OF  THE  MATRIX  PASSED  TO  THE  SSVDC  ROUTINE. 

V  :  ORTHOGONAL  MATRIX  WHOSE  COLUMNS  ARE  THE  RIGHT  SINGULAR 

VECTORS  OF  THE  MATRIX  PASSED  TO  THE  SSVDC  ROUTINE. 

WORK  :  VECTOR  USED  IN  THE  SS\’^DC  ROUTINE. 


•  INITIALIZE  THE  COUNTING  VARIABLES 

NOZERO  =  0 

NIDE  =  0 

NM  =  100 

NZB  =0 

STEP  =  1 

• 

•  FREQUENTLY  USED  FORMAT  STATEMENTS 

10  FORMAT  (A) 

50  FORMAT  (100(F12.5)) 

• 

WRITE  (6.100) 

100  FORMAT 

/'••  MODIFIED  OR  CHAINED  AGGREGATION 

. . . . . . . 


INPUT  THE  SYSTEM  MATRICES  • 


NAME  =  '"C"  MATRIX.’ 

CALL  INPUT(C.NRC,NCC.NAME.CFRMT.FCNAME) 

•  DETERMINE  THE  SMALLEST  OF  THE  TWO  DIMENSIONS 
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NMIN  =  MIN(NRC.NCC) 

* 

•  INITIALIZE  IHE  SIZE  OF  THE  SYSTEM  AND  THE  SIZE  OF  A 
SOS  =  NCC 

SIZA  =  SOS 


NAME  =  '"A"  MATRIX.’ 

CALL  INPDT(A. SOS. SOS. NAME, AFRMT.FANAME) 

• 

*  SELECT  THE  NORM  TO  BE  USED  FOR  THE  A  MATRIX. 

625  WRITE  (6.650) 

650  FORMAT  (/'ENTER  THE  NORM  OF  THE  SYSTEM  MATRIX.  A.  ' 
r  1 .  LARGEST  ELEMENT  OF  A  ' 

/’  2.  TWO-NORM  (LARGEST  SINGULAR  VALUE)' 

/’  3.  FROBENIUS  NORM' 

/ ’  ENTER  1 ,  2 .  OR  3  > ' . i) 


READ 

(5, 

•) 

I 

IF 

(I 

.EQ. 

3) 

GO 

TO 

800 

IF 

(I 

.EQ. 

2) 

GO 

TO 

750 

IF 

(I 

•  NE, 

1) 

GO 

TO 

625 

•  FIND  THE  LARGEST  ELEMENT  OF  A  AND  USE  THIS  AS  THE  NORM. 

* 

NORM  =  0.0 

DO  700  J=l.  SOS 
DO  700  1=1,  SOS 

IF  (ABS(A(I.J))  .GE.  NORM)  NORM  =  ABS(A(I.J)) 
700  CONTINUE 
GO  TO  850 


•  COMPUTE  THE  TWO  NORM 


CALL 

DUP(A.DMB,SOS.SOS) 

750 

CALL 

SSVDC(DMB . NM, SOS , SOS , SIGMA, E, D, NM. V, NM, WORX, 00 . lERR) 

IF 

(lERR  .EQ.  0)  GO  TO  760 

WRITE 

(6.755) 

755 

FORMAT 

(/'SINGULAR  VALDES  OF  A  NOT  COMPUTED  CORRECTLY' 

/'SEE  LINPACK  MANUAL  FOR  CASE  WHEN  (lERR  .NE.  0)’) 

STOP 

760 

NORM 

=  SIGMA(l) 

% 

GO  TO 

850 

• 

• 

COMPOTE  THE  FROBENIUS  NORM 

800 

NORM 

=  0 

DO  825  1=1.  SOS 

DO  825  J=l,  SOS 

NORM  =  NORM  +  A(I.J)  •Ad.  J) 
825  CONTINUE 
NORM 


I 


SQRT(NORM) 


850  WRITE  (6.860)  NORM 

860  FORMAT  (/'THE  NORM  OF  A  EQUALS:  ',F12.3) 

* 

•  CALCULATE  SIZCHK:  ELEMENTS  SMALLER  THAN  SIZCHK 

*  ARE  TREATED  AS  ZEROES.  (MACHINE  EPSILON  =  2.0**(-24)) 

« 

SIZCHK  -  ABS(NORM*((2.0)**{-24))) 

NAME  =  '"B"  MATRIX.' 

CALL  INPUT (B.NRB.NCB. NAME, BFRMT.FBNAME) 

900  CONTINUE 


*  CHOOSE  BETWEEN  PERFORMING  MODIFIED  OR  CHAINED  * 


FCHND  =  0 
950  WRITE  (6.1000) 

1000  FORMAT  (/'DO  YOU  WANT  TO  PERFORM  MODIFIED' 
/'OR  CHAINED  AGGREGATION?' 

/'  1)  MODIFIED' 

/'  2)  CHAINED' 

/'ENTER  1  OR  2  >'.i) 

« 


READ 

(5,10)  ANS 

IF 

(ANS  .EQ.  '1') 

GO  TO 

2200 

0 

IF 

(ANS  .EQ.  '2') 

GO  TO 

2000 

WRITE 

(6,1100) 

1100 

FORMAT 

(/'INCORRECT  ENTRY.  '. 

//) 

• 

GO  TO  950 

2000 

CONTINUE 

•  SET  THE  CHAINED  FLAG 

FCHND  =  1 

•  STORE  THE  NUMBER  OF  COLUMNS  IN  THE  ORIGINAL  B 

•  IN  A  TEMPORARY  NAME 

NCBTMP  =  NCB 
NCB  =1 

• 

•  CREATE  THE  NEEDED  NULL  B  MATRIX  FOR  CHAINED 
DO  2100  1=1. NRB 

B(I,1)  =  0.0 
2100  CONTINUE 
2200  CONTINUE 


•  COMPRESS  THE  COLUMNS  OF  THE  C  MATRIX  • 


•  EVENTUALLY  WE  WILL  ONLY  NEED  TO  HAVE  THE  V  MATRIX  RETURNED. 

•  FOR  NOW  WE  WANT  TO  HAVE  U  AND  V  AVAILABLE  FOR  CHECKING  PURPOSES. 


CALL  SSVDC(C.NM.NRC.NCC.SIGMA,E.D.NM,V,NM.WOR£,ll ,IERR) 
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COMPUTE  TBE  NUMBER  OF  NON-ZERO  SINGULAR  VALUES. 
DISPLAY  THE  U  AND  V  MATRIX  IF  YOU  WISH. 

NAME  =  'THE  C  MATRIX' 

CALL  USV(U. SIGMA. V.SIZCHK.NMIN.KZ.NRC.NCC, NAME) 

CHECK  THE  SIZE  OF  NZ  FOR  THE  TRIVIAL  CASE. 

IF  (NZ  .EQ.  0)  GO  TO  60600 


INITIALIZE  THE  TRANSFORMATION  MATRIX  • 


DO  4000  1=1. SOS 

DO  4000  J=l.SOS 
T(I.J)=V(J.I) 

4000  CONTINUE 

WRITE  (6.4200) 

4200  FORMAT  (/'WOULD  YOU  LIKE  TO  SEE  THE  T  MATRIX? '.i) 

READ  (5.10)  ANS 

IF  ((ANS  .NE.  'Y'). AND. (ANS  .NE.  '7’))  GO  TO  4500 
CALL  OUTPUT(T.SOS.SOS.SIZCHK) 

4500  CONTINUE 

•  CHECK  IF  THE  "C"  MATRIX  HAS  RANK  KtUAL  TO  THE  ORDER  OF  THE  SYSTEM. 

0 

IF  (NZ  .EQ.  SOS)  GO  TO  60700 


CALCULATE  THE  A12  SUBMATRIX  • 


•  SET  THE  A12  FLAG  EQUAL  TO  ZERO 

FA12  =  0 
DO  5000  1=1. NZ 

DO  5000  J=1.SIZA 
TA12(I.J)  =  0.0 

DO  5000  K=1.SIZA 

TA12(I.J)  =  TA12(I.J)  +  V(K.I)  •  A(K.J) 

5000  CONTINUE 

DO  5500  1=1. fC 

DO  5500  J=1.SIZA  -  NZ 
A12(I.J)  =  0.0 

JD  =  NZ  +  J 
DO  5500  K=1.SIZA 

A12(I.J)  =  A12(I.J)  +  TA12(I.K)  •  V(K,JD) 

IF  (ABS(A12(I.J))  .GE.  SIZCHE)  FA12  =  1 

5500  CONTINUE 

• 

WRITE  (6.5600) 

5600  FORMAT  (/'WOULD  YOU  LIKE  TO  SEE  THE  A12  SUBMATRIX?  >'.i) 
READ  (5.10)  ANS 


IF  ((ANS  .NE.  'Y').AND.(ANS  .NE.  'y'))  GO  TO  5850 
CALL  ODTPUT(A12.NZ.SIZA-NZ.SIZCHK) 

• 

5850  CONTINUE 
•  WAS  A12  =  0  (FA12  =  0) 

IF  (FA12  .EQ.  0)  GO  TO  50000 

5900  CONTINUE 

« 


•  COMPOTE  THE  NEW  B  MATRIX  • 


* 

DO  6000  I=1.SIZA 
DO  6000  J=1.NCB 
Gd.J)  =  0.0 
DO  6000  K=1.SIZA 

Gd.J)  =  Gd.J)  +  V(K.I)  •  B(NIDE+K.J) 

6000  CONTINUE 

IF  (FCHND  .EQ.  1)  GO  TO  6450 

• 

WRITE  (6.6050) 

6050  FORMAT  (/'WOULD  YOU  LIKE  TO  SEE  THE  NEW  PORTION  OF  B?  >',i) 
READ  (5.10)  ANS 

IF  ((ANS  .NE.  'Y'). AND. (ANS  .NE.  'y'))  GO  TO  6150 
CALL  OUTPDT(G.SIZA.NCB.SIZCHK) 

• 

6150  CONTINUE 

•  STORE  THIS  UPDATED  SECTION  OF  THE  B  MATRIX 

•  IN  THE  APPROPRIATE  LOCATION  OF  THE  ENTIRE  B  MATRIX. 

• 

DO  6200  I=1.SIZA 
ID  =  NIDE  +  I 
DO  6200  J=1.NCB 
BdD.J)  =  Gd.J) 

6200  CONTINUE 
• 

WRITE  (6,6300) 

6300  FORMAT  (/'WOULD  YOU  LIKE  TO  SEE  THE  ENTIRE  B  MATRIX?  > ’ , t) 
READ  (5,10)  ANS 

IF  ((ANS  .NE.  'Y') .AND.(ANS  .NE.  'y'))  GO  TO  6450 
CALL  OUTPUT(B.SOS.NCB.SIZCHK) 

• 

6450  CONTINUE 

•  CREATE  A  DUMMY  MATRIX  COMPOSED  OF 

•  THE  FIRST  NZB+NZ  ROWS  OF  B,  STARTING 

•  WITH  THE  FIRST  NON-ZERO  ROW  (NOZERO+1) . 

• 

DO  6500  1=1. NZB+NZ 

ID  =  NOZERO  +  I 
DO  6500  J=1,NCB 

GDMSd.J)  =  BdD.J) 


6500  CONTINUE 

IF  (FCHND  .EQ.  1)  GO  TO  6750 

• 

WRITE  (6.6600) 

6600  FORMAT  (/'WOULD  YOU  LIKE  TO  SEE  IHE  SUBMATRIX' 

/'OF  B  WHICH  IS  TO  BE  DECOMPOSED?  > ' . i) 

READ  (5,10)  ANS 

IF  ((ANS  .NE.  'Y'). AND. (ANS  ,NE.  'y'))  GO  TO  6750 
CALL  ODTPDT(GDMB,NZB+NZ.NCB,SIZCHK) 

* 

6750  CONTINUE 

NRGDMB  =  NZB  +  tO. 

NMINB  =  MIN ( NRGDMB. NCB) 

* 

CALL  DUP(GDMB.DMB. NRGDMB. NCB) 

•  PASS  THE  DMB  TO  SSVDC  TO  COMPRESS  ITS  ROWS. 

•  ONLY  RETURN  THE  U  MATRIX.  WE  DO  NOT  WANT  TO  OVER- 

•  WRITE  THE  CURRENT  V  MATRIX. 

CALL  SSVDC(DMB.NM, NRGDMB, NCB. SIGMA, E.U.NM.V.NM, WORK. 10, lERR) 

IF  (FCHND  .EO.  0)  GO  TO  6800 

NZB  =0 

•  SKIP  B  AND  T  UPDATE 

GO  TO  10465 

• 

6800  CONTINUE 

•  COMPUTE  THE  NUMBER  OF  NON-ZERO  SINGULAR  VALUES. 

•  DISPLAY  THE  U  MATRIX  IF  DESIRED. 

NAME  =  'A  B  SUBMATRIX. ' 

« 

CALL  USV ( U, SIGMA. V, SIZCHK, NMINB . NZB, NRGDMB , NCB , NAME) 

•  CHECK  TO  SEE  IF  GDMB  WAS  EQUAL  TO  ZERO. 

•  IF  SO,  THEN  SKIP  THE  B  AND  T  UPDATE. 

IF  (NZB  .EQ.  0)  GO  TO  10465 


UPDATE  THE  B  MATRIX  • 


*  ZERO  OUT  THE  APPROPRIATE  ELEMENTS  IN  THE  B  MATRIX. 

*  THIS  SUBMATRIX  IS  THE  SIZE  OF  GDMB. 

DO  10000  1=1. NRGDMB 
ID  =  NOZERO  +  I 
DO  10000  J=1,NCB 
B(ID.J)  =  0.0 
10000  CONTINUE 


^  i 


S7 


• 

•  PERFORM  THE  TRANSFORMATION  ON  THE  B  MATRIX 

•  WITH  IHE  D  MATRIX,  COMPRESSING  THE  ROWS  OF 

•  THE  NEW  B  DOWN. 

DO  10100  1*1, NZB 

ID  *  NOZERO  +  NRGDMB  -  NZB  +  I 
DO  10100  J=1,NCB 

DO  10100  E=l, NRGDMB 

B(ID,J)  =  B(ID,J)  +  D(K,I)  •  GDMB(K,J) 

10100  CONTINUE 

WRITE  (6,10150) 

10150  FORMAT  (/’WOULD  YOU  LIKE  TO  SEE  THE  TRANSFORMED  B  MATRIX' 

/•AFTER  THE  APPROPRIATE  ROWS  HAVE  BEEN  COMPRESSED?  >',t) 
READ  (5,10)  ANS 

IF  ((ANS  .NE.  'Y'). AND. (ANS  .NE.  'y'))  GO  TO  10165 
CALL  OUTPDT(B.SOS.NCB,SI2CHK) 

10165  CONTINUE 


UPDATE  THE  TRANSFORMATION  MATRIX  • 


•  THIS  MUST  BE  DONE  IN  TWO  STAGES.  BECAUSE  WE 

•  COMPRESS  THE  ROWS  DOWN  IN  THE  B  SUBMATRIX. 

« 

•  CALCULATE  THE  UPPER  PART  OF  THE  NEW  T  MATRIX. 

DO  10210  I=1,NRGDMB-NZB 
ID  =  NZB  +  I 
DO  10210  J=l.SOS 
TDMB(I.J)  =  0.0 
DO  10210  K=l. NRGDMB 

TDMB(I.J)  =  TDMBd.J)  +  U(K.ID)  •  T(NOZERO+K, J) 

10210  CONTINUE 

• 

•  CALCULATE  THE  LOWER  PART  OF  IHE  NEW  T  MATRIX. 

• 

DO  10220  1=1, NZB 

ID  =  NRGDMB  -  NZB  +  I 
DO  10220  J»1 . SOS 

TDMB(ID,J)  =  0.0 
DO  10220  K=l, NRGDMB 

TDMBdD.J)  =  TDMBdD.J)  +  D(K,I)  •  T(K+NOZERO. J) 

10220  CONTINUE 

WRITE  (6.10250) 

10250  FORMAT  (/'WOULD  YOU  LIKE  TO  SEE  THE  UPDATED  ROWS  IN  T' 
/’DUE  TO  THE  D  GENERATED  WHILE  COMPRESSING  THE’ 

/ 'ROWS  OF  B7  >’ .1) 

READ  (5,10)  ANS 

IF  ((ANS  .NE.  'Y') .AND. (ANS  .NE.  'y'))  GO  TO  10265 
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CALL  OUTPDT(TDMB.NRGDMB.SOS.SIZCHK) 

10265  CONTINUE 

« 

•  STORE  TOMB  IN  THE  APPROPRIATE  LOCATION  WITHIN  T 

« 

DO  10300  I=1,NRGDMB 
ID  =  NOZERO  +  I 
DO  10300  J=l.SOS 
T(ID.J)  =  TOMB(I.J) 

10300  CONTINUE 

WRITE  (6.10350) 

10350  FORMAT  (/'WOULD  YOU  LIKE  TO  SEE  THE  NEW  T?  >',i) 

READ  (5,10)  ANS 

IF  ((ANS  .NE,  'Y') .AND. (ANS  .NE,  'y'))  GO  TO  10365 
CALL  OCTPUT(T.  SOS.SOS.SIZCHK) 

* 

10365  CONTINUE 


UPDATE  TOE  A12  SUBMATRIX  WITH  THIS  U  • 
STATE  SPACE  TRANSFORMATION.  • 


*  A12  WILL  BE  DIVIDED  INTO  TWO  PARTS,  A12B  AND  A12S. 

*  CALCULATE  A12B. 

DO  10400  1=1, NZB 

DO  10400  J=1,SIZA-NZ 
A12B(I.J)  =  0.0 
DO  10400  K=1,NRGDMB 

A12B(I,J)  =  A12B(I,J)  +  D(K.I)  •  A12(K,J) 

10400  CONTINUE 

WRITE  (6,10450) 

10450  FORMAT  (/'WOULD  YOU  LIKE  TO  SEE  TOE  A12B  SUBMATRIX?  >',i) 
READ  (5,10)  ANS 

IF  ((ANS  .NE.  'Y') .AND. (ANS  .NE.  'y'))  GO  TO  10465 
CALL  OUTPUT(A12B.NZB,SIZA-NZ.SIZCHK) 

* 

10465  CONTINUE 

* 

*  CALCULATE  A12S. 

•  SET  TOE  A12S=0  FLAG  EQUAL  TO  ZERO. 

FA12S  =  0 

DO  10500  I=1.NEGDMB  -  NZB 
ID  =  NZB  +  I 
DO  10500  J=1,SIZA-NZ 
A12S(I,J)  =  0.0 
DO  10500  K=1,NRGDMB 

A12S(I,J)  =  A12S(I,J)  +  D(K.ID)  *  A12(K,J) 

IF  (ABS(A12S(I,J))  .GT.  SIZCHK)  FA12S  =  1 


10500  CONTINUE 


WRITE  (6,10550) 

10550  FORMAT  (/'WOULD  YOU  LIKE  TO  SEE  TBE  A12S  SUBMATRIX?  >'.i) 
READ  (5,10)  ANS 

IF  ((ANS  .NE.  'Y'). AND. (ANS  ,NE.  'y'))  GO  TO  10565 
CALL  0UrPUT(A12S,NEGDMB-NZB,SIZA-NZ,SIZCHK) 

• 

10565  CONTINUE 

IF  (FA12S  .EQ.  0)  GO  TO  50200 


•  CALCULATE  THE  A22  SUBMATRIX.  • 


•  TA12  IS  USED  AS  IHE  INTERMEDIATE  STORAGE  LOCATION, 

DO  10600  I=1,SIZA-NZ 
ID  =  NZ  +  I 
DO  10600  J=1.SIZA 
TA12(I,J)  =  0.0 
DO  10600  K=1,SIZA 

TA12(I.J)  =  TA12(I.J)  +  V(K.ID)  •  A(K,J) 

10600  CONTINUE 

• 

•  STORE  TEE  A22  SUBMATRIX  IN  THE  UPPER  LEFT  HAND 

•  CORNER  OF  THE  A  MAIRII. 

• 

DO  10700  r^l,SIZA-NZ 
DO  10700  J=1,SIZA-NZ 
JD  =  NZ  +  J 
A(I,J)  =  0.0 
DO  10700  K=1,SIZA 

A(I.J)  =  A(I.J)  +  TA12(I.K)  •  V(K.JD) 

10700  CONTINUE 

• 

WRITE  (6,10750) 

10750  FORMAT  (/’WOULD  YOU  LIKE  TO  SEE  THE  A22  SUBMATRIX?  >’,t) 
READ  (5,10)  ANS 

IF  ((ANS  .NE  'Y') .AND, (ANS  .NE.  'y'))  GO  TO  10760 
CALL  OUTPDT(A.SIZA-NZ,SIZA-NZ,SIZCHK) 

10760  CONTINUE 

• 

•  UPDATE  THE  SIZA  VARIABLE  TO  BE  THE  SIZE  OF  *HE 

•  NEW  A22  SUBMATRIX. 

• 

SEZA  =  SIZA  -  NZ 

•  UPDATE  THE  VARIABLE  REPRESENTING  THE  NUMBER  OF 

•  IDENTITY  ELEMENTS.  (THIS  VARIABLE  IS  USED  DURING 

•  THE  T  UPDATE  WITH  V.  ) 


NIDE  =  NIDE  +  NZ 
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IDEM  I!  ICATION  OF  STEP  TO  DSER  • 


•  INDICATE  TO  DSER  WHAT  STEP  IN  THE  DECOMPOSITION 


♦  HE 

IS  IN. 

NAME 

B  *  * 

IF 

(FCHND  .EQ.  0)  NAME  >  'MODIFIED' 

WRITE 

(6,10770)  STEP  .  NAME 

10770 

FORMAT 

/'THIS  COMPLETES  '.13.'  STAGE(S)  OF  '.AS, 
/'CHAINED  AGGREGATION.') 

WRITE 

(6,10772)  NIDE  ,  NIDE 

10772 

FORMAT 

(/'THE  AGGREGATE  SYSTEM  IS  NOW  '.13,'  z  ',13. 

IF 

(STEP  .NE.  1)  GO  TO  10780 

WRITE 

(6.10774) 

10774 

FORMAT 

(/'WOULD  YOU  LIRE  THIS  INFORMATION  TO  BE' 
/'WRITTEN  OUT  TO  A  FILE?  >'.t) 

READ 

(5,10)  ANS 

IF 

((ANS  .EQ.  'Y'). OR. (ANS  .EQ.  'y'))  GO  TO  10775 

INFO 

=  .FALSE. 

GO  TO  10780 

10775 

INFO 

=  .TRUE. 

WRITE 

(6,10776) 

10776 

FORMAT 

(/'ENTER  THE  FILE  NAME  >',i) 

READ 

(5,10)  INAME 

OPEN 

(UNIT=3 ,FILE=INAME) 

REWIND 

3 

WRITE 

(3,10778) 

10778 

FORMAT 

(/'THIS  FILE  CONTAINS  THE  INFORMATION  ON' 

/'HOW  THE  AGGREGATE  SYSTEM  GROWS  WITH' 

/'EACH  STEP  OF  MODIFIED  CHAINED  AGGREGATION.', 

CLOSE 

(UNIT=3) 

10780 

IF 

(.NOT.  INFO)  GO  TO  10790 

OPEN 

(DNIT=3 , FILE* INAME) 

WRITE 

(3,10785)  STEP  ,  NIDE  ,  NIDE 

10785 

FORMAT 

(/'AFTER  ',13,'  STAGE(S)  OF  MODIFIED  CHAINED' 

/'AGGREGATION.  THE  AGGREGATE  SYSTEM  IS' 
/,I3,'  X  ',13) 


10790  CONTINDE 

STEP  *  STEP  +  1 


COMPRESS  THE  COLUMNS  OF  A12S  * 


WILL  RETURN  BOTH  U  AND  V  MATRICES  FOR  NOW. 
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« 


NBA12S  -  NBGDMB  -  ICB 
*  RECALL  NCA12S  -  SIZA. 

CALL  SSVDC(A12S.NM.NRA12S. SIZA. SIGMA. E.D.NM.V.NM, WORK, 11, lERR) 


NAME  >  'THE  A12S  SUBMAISIX. ' 

NMIN  -  MIN(NRA12S,SIZA) 

CALL  IJSV(U.  SIGMA.V,  SIZCHK, NMIN, fC.NRAllS,  SIZA,  NAME) 


*  IF  NZ  EQUALS  SIZA  THEN  THE  SYSTEM  WILL  NOT  AGGREGATE. 

*  THIS  IS  BECAUSE  A12S  HAS  FULL  COLUMN  RANK. 

IF  (NZ  .EQ.  SIZA)  GO  TO  60800 


•  UPDATE  T  • 


•  DO  THE  UPDATE  IN  TWO  STEPS. 

•  COMPUTE  AFFECTED  ROWS  OF  T. 

• 

DO  10800  I>:1 ,  SIZA 

DO  10800  J^l.SOS 
TDMBd.J)  «  0.0 
DO  10800  E>°1 .  SIZA 

TDMBd.J)  =  TDMBd.J)  +  V(K,I)  •  T(NIDE+K.J) 

10800  CONTINUE 

* 

•  STORE  THESE  AFFECTED  ROWS  IN  THE  APPROPRIATE  LOCATION 

•  WITHIN  T. 

• 

DO  10900  I-l . SIZA 

ID  >  NIDE  I 
DO  10900  J-l.SOS 
TdD.J)  -  TDMBd.J) 

10900  CONTINUE 
* 

WRITE  (6,10950) 

10950  FORMAT  (/'WOULD  YOU  LIKE  TO  SEE  THIS  NEW  T* 

/'AFTER  BEING  UPDATED  WITH  V7  >',i) 

READ  (5,10)  ANS 

IF  ((ANS  .NE.  'V') .AND. (ANS  .NE.  'y'))  GO  TO  10970 
CALL  (XrrPUT(T,SOS,SOS.SIZCHK) 

10970  CONTINUE 


COMPUTE  THE  NEW  A12  SOBMATHIl  * 


SET  THE  A12-0  FLAG  TO  ZERO 
FA12  -  0 
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•  CALCULATE  IHE  UPPER  HALF  OF  IHE  A12  SUBNATRIX. 

• 

DO  11000  I-1,)CB 

DO  11000  J>1,SIZA-NZ 
JD  =  NZ  +  J 
A12(I.J)  =  0.0 
DO  11000  .  SIZA 

A12(I,J)  -  A12(I.J)  +  A12B(I,K)  •  V(K,JD) 

IF  (ABS(A12(I.J))  .GT.  SIZCHX)  FA12  »  1 

11000  CONTINUE 

* 

•  CALCULATE  THE  LOWER  HALF  OF  THE  A12  SUBMATRIX  IN 

•  TWO  SECTIONS. 

« 

DO  11100  1=1.  IC 

DO  11100  J=1.SIZA 
TA12(I.J)  =  0.0 
DO  11100  C=1 . SIZA 

TA12(I.J)  =  TA12(I,J)  +  V(K,I)  *  A(K.J) 

11100  CONTINUE 

DO  11200  I-1,^E 

ID  =  NZB  +  I 

DO  11200  J=1,SIZA  -  NZ 
JD  =  NZ  +  J 
A12<ID,J)  =  0.0 
00  11200  E=1.SIZA 

A12(ID.J)  -  A12(ID,J)  +  TA12(I,K)  •  V(K,JD) 

IF  (ABS(A12(ID.J))  .GT.  SIZCHE)  FA12  =  1 

11200  CONTINDE 

« 

WRITE  (6.11250) 

11250  FORMAT  (/'WOULD  YOU  LIEE  TO  SEE  THE  NEW  A1 2  SUBMATRIX?  >'.i) 
READ  (5,10)  ANS 

IF  ((ANS  .NE.  'Y'). AND. (ANS  .NE.  'y'))  GO  TO  11270 
CALL  0DTPDT(A12,ICB+NZ,SIZA-NZ. SIZCHE) 

11270  CONTINUE 

« 

•  CHECE  TO  SEE  IF  A12=0.  IF  SO.  SYSTEM  AGGREGATES. 

IF  (FA12  .EQ.  0)  GO  TO  50000 


UPDATE  THE  VARIABLE  REPRESENTING  TEE  NUMBER  OF  ZERO 
ROWS  IN  THE  B  MATRIX. 

NOZERO  =  NIDE  -  NZB 


CONTINUE  TEE  PROCESS  OF  MODIFIED  CHAINED  AGGREGATION.  • 


1 


40900  GO  TO  5900 
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PROGRAM  EXITS  • 


50000  WRITE  («, 50100) 

50100  FORMAT  (/’PROGRAM  EXIT.  SYSTEM  AGGREGATES.' 
/’THE  A12  SUBMATRIX  -  0.’) 

GO  TO  70000 

50200  WRITE  (6,50250) 

50250  FORMAT  (/’PROGRAM  EXIT.  SYSTEM  WILL  AGGREGATE. 
/’THE  MATRIX  A12S  -  0.’) 

GO  TO  70000 


PROGRAM  EXITS  ASSOCIATED  WITH  NO  AGGREGATION.  • 


60600  WRITE  (6,60650) 

60650  FORMAT  (/’HtOGRAM  EXIT.  TRIVIAL  CASE’ 

/’THE  C  MATRIX  IS  ZERO.’) 

GO  TO  80000 

* 

60700  WRITE  (6,60750) 

60750  FORMAT  (/’PROGRAM  EXIT.  SYSTEM  WILL  NOT  AGGREGATE. ’ 

/’THE  C  MATRIX  HAS  A  RANK  EQUAL  TO  THE’ 

/’DIMENSION  OF  THE  SYSTEM.’) 

GO  TO  70000 

« 

60800  WRITE  (6,60850) 

60850  FORMAT  (/'PROGRAM  EXIT.  SYSTEM  WILL  NOT  AGGREGATE. ’ 

/’THE  A12S  SUBMATRIX  HAS  FULL  COLUMN  RANK.’) 

GO  TO  70000 

• 

70000  CONTINUE 

WRITE  (6,70100) 

70100  FORMAT  (/’WOULD  YOU  LIKE  TO  SEE  THE  TRANSFORMED' 

/’SYSTEM  MATRICES?  >’,i) 

READ  (5,10)  ANS 

IF  ((ANS  .NE.  ’Y’) .AND. (ANS  .NE.  'y’))  GO  TO  80000 

*  COMPUTE  THE  TRANSFORMED  A  MATRIX. 

CALL  NULT(T,0,AFBMT,FANAME, SOS, SOS, A.1,AFRMT,FANAME, SOS, SOS, TOMB) 
DO  70200  I-l,SOS 
DO  70200  J-l,SOS 
DMB(I,J)  -  T(J.I) 

70200  CONTINDE 

• 

CALL  NDLT(TDMB,0,AFRMT,FANAME,SOS,SOS,DMB,0,AFRMT,FANAME,SOS, 

SOS, A) 
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WRITE  (6,70300) 

70300  FORMAT  (/'THIS  IS  THE  TRANSFORMED  "A"  MATRIX:') 

CALL  0(ITPUT(A.S0S,S0S,SIZCHK) 

•  COMPUTE  THE  TRANSFORMED  "B”  MATRIX. 

* 

IF  (FCHND  .EQ.  1)  NCB  NCBTMP 

CALL  MDLT(T.  0  .BFRMT,  FBNAME,  SOS,  SOS,  B,  1  ,BFRMT,FBNAME,  SOS,  NCB ,  G) 
WRITE  (6,70400) 

70400  FORMAT  (/'THIS  IS  THE  TRANSFORMED  "B"  MATRIX.') 

CALL  (M7TPDT(G,SOS,NCB,SIZCHE) 

•  COMPUTE  THE  TRANSFORMED  "C"  MATRIX. 

•  (RECALL  DMB  IS  THE  TRANSPOSE  OF  THE  "T"  MATRIX.) 

CALL  MULT(C,1,CFRMT,FCNAME,NRC,NCC,DMB,0,CFRMT,FCNANE, 
SOS,SOS,TDMB) 

WRITE  (6,70500) 

70S00  FORMAT  (/'THIS  IS  THE  TRANSFORMED  "C  MATRIX.’) 

CALL  OOTPUT(TDMB,NRC,NCC,SIZCHX) 

• 

WRITE  (6,70600) 

70600  FORMAT  (/'WOULD  TOD  LIIE  TO  SEE  THE  FINAL  "T”  MATRIX?  >’,i) 

READ  (3,10)  ANS 

IF  ((ANS  .NE.  'T') .AND. (ANS  .NE.  'y'))  GO  TO  70800 

• 

CALL  OUTPUT(T,SOS.SOS.SIZCHX) 

70800  WRITE  (6,70900) 

70900  FORMAT  (/'WOULD  YOU  LIKE  TO  SEE  THE  FINAL  "T"  TRANSPOSE' 
/'MATRIX?  >',l) 

READ  (3,10)  ANS 

IF  ((ANS  .NE.  'T') .AND. (ANS  .NE.  'y'))  GO  TO  80000 

CALL  OUTPUT(DMB.SOS,SOS.SIZCBK) 

80000  CONTINUE 
STOP 
END 


•  SUBROUTINES  USED  IN  PROGRAM  ABOVE  • 


« 


•  SUBROUTINE  USV  • 


•  THIS  SUBROUTINE  COMPUTES  THE  NUMBER  OF  NON-ZERO 


y 


§5 


•  SINGULAR  VALDBS  ASSOCIATED  WITH  A  DECOMPOSED 

•  MATRIX.  IT  ALSO  ASKS  IP  TOU  WOULD  LUE  TO  SEE 

•  THE  U  AND/OR  V  MATRIX  CREATED  DURING  TEE  DECOMPOSITION. 

•  DEPENDING  ON  WHERE  IN  THE  PROGRAM  USV  IS  CALLED, 

• 

SUBROUTINE  USV(U. SIOMA.V.SIZCHK.MIN.N.NR.NC.TIFE) 
CHARACIER*!  ANS 


CHARACTER*20 

TYPE 

REAL 

U(IOO.IOO) 

.  V(IOO.IOO) 

,  SIGMAdOO) 

REAL 

SIZCHK 

INTEGER 

MIN 

.  N 

.  NR 

INTEGER 

I 

.  lERR 

•  FOEMAT  STAIEMENTS  USED 

910 

FORMAT 

(A) 

950 

% 

FORMAT 

(100(F12.5)) 

WRITE 

(6.990)  TYPE 

990 

FORMAT 

(/'THIS  DECOMPOSITION  IS  OF 

NC 


WRITE  (6,1000) 

1000  FORMAT  (/'WOULD  TOU  LIKE  TO  SEE  THE  U  MATRIX?  >',t) 
READ  (5,910)  ANS 

IF  ((ANS  .NE.  'T'). AND. (ANS  .NE.  'y'))  GO  TO  2050 

• 

WRITE  (6.1500) 

1500  FORMAT  (/'THE  OORBESPOND1N6  U  MATRIX  IS:') 

CALL  OSFLAX(U.NR.NR) 


2050 

« 

2100 


2200 

• 

2300 

2350 


2400 


CONTINUE 
WRITE  (6,2100) 

FORMAT  (/'THE  SINGULAR  VALUES  ARE:') 

WRITE  (6,950)  (SIOMA(I) .I«1,MIN) 

CHECK  THE  MAGNITUDE  OF  THE  SINGULAR  VALUES 
N^O 

DO  2200  I-1,MIN 

IF  (SIOMA(I).LT.SIZCHK)  GO  TO  2300 

I^Nfl 

CONTINUE 

WRITE  (6,2350)  N 

FORMAT  (/'THE  NUMBER  W  INDEPENDENT  COLUMNS  IN' 

/'THE  mCOMPOSED  MATRIX  IS:  ',15) 

IF  (TYPE  .EQ.  'A  B  SUBMATRIX.')  GO  TD  3050 
WRITE  (6,2400) 

FORMAT  (/'WOULD  YOU  LIEE  TO  SEE  THE  V  MATRIX?  >',!) 
READ  (5,910)  ANS 

IF  ((ANS  .NE.  'Y'). AND. (ANS  .NE.  'y'))  GO  TO  3050 


WRITE  (6,2500) 


2S00 

FORMAT 

(/'THE  CORRESPONDING  V  MATRIX  IS:') 

CALL 

DSPLAY(V.NC.NC) 

SOSO 

CONTINUE 

WRITE 

(6,3500)  lERR 

3500 

• 

FORMAT 

RETURN 

END 

(/'lERR  EQUALS  ',14) 

SUBROUTINE  DUP  • 


THIS  SUBROUTINE  DUPLICATES  A  MATRIX.  THIS  IS  NEEDED 
SINCE  THEN  A  MATRIX  IS  PASSED  TO  SSVDC  IT  RETURNS  IN 
AN  ALTERED  FORM. 


SUBROUTINE  DUP( SOURCE. DUMMY, NROW.NCOL) 


• 

REAL 

SOURCEdOO.lOO) 

.  DUMMYdOO.lOO) 

INTEGER 

NROW 

,  NCOL 

INTEGER 

* 

I 

,  J 

DO  1000  I>=l.NROW 

DO  1000 

J-l,NCOL 

DUMMYd.J)  »  SOURCEd.J) 

1000  CONTINUE 
RETURN 
END 


SUBROUTINE  OUTPUT  • 


THIS  SUBROUTINE  ALLOWS  THE  USER  TO  SEE  THE 
PARTICULAR  MATRIX  ON  THE  SCREEN  IN  EITHER  F6.3 
FORMAT  OR  SIMPLY  AS  X's  AND  O's  (IF  THE  STRUCTURE 
IS  ALL  THAT  IS  DESIRED) .  THE  USER  MAY  ALSO  OUTPUT 
THE  MATRIX  TO  A  FILE  IN  EITHER  OF  THESE  FORMATS. 
THE  USER  MUST  SUPPLY  THE  FILENAME  TO  BE  USED. 

SUBROUTINE  OUTPUT(X , NROW , NCOL, SIZCHK) 


CHARACTER*! 

ANS 

CHARACTER*! 

XCdOO.lOO) 

CHARA(TER*20 

FNAME 

REAL 

XdOO.lOO) 

,  siza 

INTEGER 

I 

.  J 

INTEGER 

NROW 

,  NCOL 

INTEGER 

FORM 

LOGICAL 

Four 

FREQUENTLY  USED 

FORMAT  STATEMENTS. 

A 


fT 


10 

FOBMAT 

(8(F9.3,1X)) 

15 

FOBNAT 

(7(B15.7.1X)) 

20 

FOBMAT 

(A) 

30 

0 

FOBXAT 

e 

► 

• 

H 

50 

WBITE 

(6.100) 

100 

FOBMAT 

(/'WHAT  FOBMAT  DO  TOO  WISH  TO 

SEE' 

/'THIS  MATBIX  IN?' 

/'  1)  NDMBBICAL  VALDES' 

/'  2)  STBDCIDBE  ONLT  (X  AND  0) ' 

• 

/'ENTEB  1  OB  2  >'.<) 

BRin 

(5.150)  FOBM 

150 

0 

FOBMAT 

(11) 

WBITE 

(6.200) 

200 

FOBMAT 

(/'WOULD  TOD  LIKE  THIS  MATBIX 

STOBED' 

• 

/'IN  A  FILE?  >'.*) 

BRAD 

(5.20)  ANS 

IF 

((ANS  .SO.  'T'). OB. (AMS  .EO.  ‘ 

y'))  G4 

FOOT  » 

.FALSE. 

GO  TO  400 

275 

WBITE 

(6.300) 

300 

FOBMAT 

(/'ENTEB  THE  FILENAME  FOB  THIS  MATBIX 

BEAD 

(5.20)  FNAME 

FOOT  - 

.TBDE. 

400 

IF 

(FOBM  .BQ.  1)  GO  TO  500 

IF 

(FOBM  .EQ.  2)  GO  TO  1000 

WBITE 

(6.450) 

450 

FOBMAT  (/'IN(X)BBECT  ENTBT(S) . ' .//} 

GO  TO  50 

500 

CALL 

DSFLAKX,  NBOW.  NOOL) 

IF 

(.NOT.  FOOT)  BBTUBN 

600 

OPEN 

(DNIT^l .FILE-FNAME) 

BEWIND  1 

WBITE 

(1,«)  NBOW  .  NOH. 

DO  650 

I-l.NBOW 

650 


WKIIB  (1.15)  (Z(I.J),J-1.NOOL) 

CONTINDE 
CLOSE  (ONIT-l) 
lETUBN 


1000  DO  1050  I-l.NSOV 

DO  1050  j-i.Nora. 

ICd.J)  -  'O' 

IF  (ABS(X(I,J))  .GT.  SIZCHE)  XCCl.J)  -  'X' 


I 


»8 


1050  CONTINUE 

• 

IF  (FOOT)  GO  TO  1300 
DO  1100  I»l.NROW 

WRITE  (6,30)  (XC(I.J) ,J=l,NCOL) 
1100  CONTINUE 
RETURN 

• 

1300  OPEN  (UNIT=1.FILE*=FNAME) 

REWIND  1 

DO  1400  I-l.NROW 

WRITE  (6,30)  (XC(I,J).J=l.NCOL) 
WRITE  (1,30)  (XC(I,J) ,J=l,NCOL) 
1400  CONTINUE 

CLOSE  (UNITED 

RETURN 

END 


SUBROUTINE  INPUT  • 


TBIS  SUBROUTINE  READS  IN  THE  APPROPRIATE  MATRIX  FROM 
A  DATA  FILE.  THE  MATRIX  MAT  BE  STORED  IN  NORMAL  FORM 
OR  AS  A  SPARSE  MATRIX.  YOU  ARE  ALSO  ALLOWED  THE  ABILITY 
TO  VERIFY  THE  MATRIX  READ  IN. 


200 

300 


250 


252 

255 


SUBROUTINE  INPDT(X,NROW,NCOL,TYPE.FRMT,FNAME) 


CHARACTER*! 

ANS 

,  FRMT 

CHARACTER*20 

FNAME 

,  TYPE 

REAL 

XdOO.lOO) 

INTEGER 

I 

.  J 

INTEGER 

NROf 

.  NCOL 

INTEGER 

NR 

.  NC 

WRITE  (6,200) 

TYPE 

FORMAT  (/'ENTER  THE  NAME  OF  THE  DATA  FILE  FOR  THE  '.All,'  > ' . i) 

READ  (5,300)  FNAME 

FORMAT  (A) 

IF  (TYPE  .EQ.  '"A"  MATRIX.')  GO  TO  252 
WRITE  (6,250)  TYPE 

FORMAT  (/'ENTER  THE  DIMENSIONS  OF  THE  ',A10.'  (ROWS  x  COLS).  >'.i) 

READ  (5.*)  NROW.NCOL 

WRITE  (6,255) 

FORMAT  (/'WHAT  FORMAT  IS  THIS  MATRIX  IN?' 

/ '  A)  SPARSE' 

/'  B)  NORMAL' 

/'  ENTER  A  OR  B  >',i) 

READ  (5.300)  FRMT  [ 

IF  ((FRMT  .EQ.  '  A' )  .OR.  (FRMT  .EQ.  '»'))  GO  TO  310  jj 

IF  ((FRMT  .EQ.  'B' ) .OR. (FRMT  .EQ.  'b'))  GO  TO  305  I 

GO  TO  252 


305  OPEN  (1INrP-l,FILE»FNAMB) 

SEVIND  1 

READ  (1.*)  NR  .  NC 

IF  ((NR  .EQ.  NSOV}.AND.(NC  .£Q.  NCOL))  60  TO  307 
WRllE  (6.306) 

306  FORMAT  (/'PROGRAM  EXIT]  MATRIX  DIMENSIONS  WHICH  WERE' 

/'INPUT  AND  THOSE  IN  THE  DATAFILE  DO  NOT  AGREE.') 

STOP 

307  READ  (!.•)  ((X(I,J)  .J-l.NCOL)  .I-l.NRW) 

CLOSE  (DNIT^l) 

GO  TO  350 


310  DO  308  I-l.NROV 

DO  308  J-l.NCOL 
X(I,J)-0.0 
308  CONTINDE 

OPEN  (DNIT-1,FILE-FNAME) 

REWIND  1 

312  READ  (1.*. END-350)  I  .  J  .  X(I,J) 

60  TO  312 

350  WRITE  (6,360)  TTPE 

360  FORMAT  (/'DO  TOD  WANT  TO  VERIFT  THE  '.AlO.'T  >',i) 
READ  (5.300)  ANS 

IF  ((ANS.NE.'I').AND.(ANS.NB. 'y'))  RBTIIRN 
WRITE  (6.400)  TtPE 

400  FORMAT  (/'THE  FOLLOWING  '.AlO.'  WAS  READ  IN:') 

CALL  DSPLATd.NROW.NOOL) 

1000  CONTINUE 

RETURN 

END 


SUBROUTINE  MOLT  * 


THIS  SUBROUTINE  MULTIPLIES  TWO  MATRICES  TOGETHER  AND 
RETURNS  THE  RESULT  IN  P.  ONE  OR  BOTH  (J¥  THE  FILES  MAI 
HAVE  TO  BE  BKAD  IN  FROM  DATAFILES,  THIS  OPTION  IS  OON- 
mJLLED  BT  A  FLAG  FOR  EACH  MATRIX. 

1  MEANS  READ  THE  FILE  IN  FROM  A  DATAFILE. 

0  MEANS  THE  MATRIX  WAS  PASSED  IN  TEE  CALL. 

•FRMT  INDICATES  THE  FORMAT  THE  PARTICULAR 
MATRIX  IS  STORED  IN. 

SUBROUTINE  MULT(X,XF,XFRMT.XNAME.NRX,NCX,7,IF.TFRMT,INAMB,NRI, 
NCT.P) 

CHAEACIER*!  XFRMT  ,  TFRMT 

CHARACTER«20  XNAME  ,  INAME 


100 


LOGICAL 

INTEGER 

INTEGER 

INTEGER 


XF 

I 

NRX 

RD 


.  TF 
.  J 
.  NRY 
,  CD 


.  E 
.  NCX 


NCT 


REAL 

X(IOO.IOO),  YdOO.lOO).  PdOO.lOO) 

IF 

(NCX  .EQ.  NRY)  GO  TO  50 

WRITE 

(6,30) 

30 

FORMAT 

(/'THE  MATRICES  ARE  NOT  COMPATABLEl ' ) 

• 

STOP 

• 

DOES  THE 

MATRIX  NEED  TO  BE  READ  IN  FROM  A  DATA  FILE? 

50 

• 

IF 

(.NOT.XF)  GO  TO  500 

WHAT  FORMAT  IS  THE  MATRIX  IN?  ( A= SPARSE. B^NORMAL) 

IF 

((XFRMT  .EQ.  'A') .OR.(XFRMT  .EQ.  'a'))  GO 

IF 

((XFRMT  .EQ.  'B') .OR. (XFRMT  .EQ.  'b'))  GO 

WRITE 

(6.100) 

100 

FORMAT 

(/'THE  X  MATRIX  IN  THE  MULT  SUBROUTINE' 

/'WAS  NOT  ASSIGNED  A  FORMAT  TYPE') 

0 

STOP 

200 

OPEN 

(UNIT^l  .FILE==XNAME) 

REWIND 

1 

READ 

(1.*)  RD  .  CD 

READ 

(1.*)  ((X{I.J).J*1.NCX).I=1,NRX) 

aosE 

(DNIT=1) 

% 

GO  TO 

500 

300 

DO  350 

1=1. NRX 

DO 

350  J=1.NCX 

350 


X(I.J) 

CONTINUE 


0,0 


OPEN  (DNIT^l,FILE*=XNAJffi) 

REWIND  1 

360  READ  (1,*.END=500)  I  ,  J  ,  X(I.J) 

60  TO  360 

•  DOES  THE  Y  MATRIX  NEED  TO  BE  READ  IN  FROM  A  DATAFILE? 

500  IF  (.NOT.YF)  GO  TO  1000 

« 

•  WHAT  FORMAT  IS  IHE  "Y"  MATRIX  IN?  (A-SPARSE.B-NORMAL) 

• 

IF  ((YFRMT  .EQ.  ' A' ) .OR. (YFRMT  .EQ.  '«'))  GO  TO  800 

IF  ((YFRMT  .EQ.  ’B') .OR. (YFRMT  .EQ.  'b'))  GO  TO  600 

WRITE  (6.550) 

550  FORMAT  (/'THE  "Y"  MATRIX  IN  THE  MULT  SUBROUTINE' 

/'WAS  NOT  ASSIGNED  A  FORMAT  TYPE') 


STOP 
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600  OPEN  (DNIT-1,FILE-1NAJ1B} 

SEffIND  1 

READ  (l.«)  RD  ,  CO 

READ  (I,*)  ((I(I,J),J-1.NCY).I-1.NRY) 

CLOSE  (DNITWI) 

GO  10  1000 

• 

800  DO  850  I-1,NRY 

DO  850  J-l.NCY 
Y(I,J)  -  0.0 
850  CONTINUE 

• 

OPEN  (DNIl^l.FILE-YNAME) 

REfflND  1 

860  READ  (1,«. END-1000)  I  .  J  ,  Y(I.J) 

GO  TO  860 

• 

1000  DO  1200  I-l.NRZ 

DO  1200  J-l.NCY 
P(I.J)  -  0.0 
DO  1200  K-l.NCX 

P(I,J)  -  P(I,J)  +  X(I,R)*Y(K.J) 

1200  CONTINUE 
RETURN 
END 


SUBROUTINE  DSFLAY  * 


•  THIS  SUBROUTINE  ALLOWS  A  MATRIX  TO  BE 

•  DISPLAYED  THE  TERMSiAL  SCREEN  IN 

•  GROUPS  OF  8  COLUMNS. 

« 

SUBROUTINE  DSPLAYCX.NROW.NOH.) 

INTEGER  I  ,  J  .  START 

INTEGER  NROr  ,  NCOL  ,  FINISH 

REAL  X( 100, 100} 

• 

START  -  1 

FINISH  -  8 

100  IF  (FINISH  .GE.  NOOL)  FINISH  -  NCOL 
WRITE  (6,150)  START  ,  FINISH 
150  FORMAT  (/'  COLS.  ',13,'  TO  ',13) 

DO  200  I-l.  NROW 

WRITE  (6,175)  (X(I,J),J-START,FINISH) 
175  FORMAT  (8(F9.3,1X)) 

200  CONTINDE 

IF  (FINISH  .GE.  NOOL)  GO  TO  300 
START  -  START  +  8 

FINISH  -  FINISH  *  8 
GO  TO  100 
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300  REITJRN 
END 


THE  FOLLOWING  RODTINES  HAVE  BEEN  TAKEN  • 
FROM  THE  LINPACK  MATHEMATICAL  SOFTWARE.  • 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SHBRODTINE  SSVDCd.LDX.N.P,  S.E.D, LDD.V.LDV, WORK,  JOB.  INFO) 

INTEGER  LDX.N.P.LDD.LDV. JOB. INFO 

REAL  X(LDX.l) .S(l) .E(l) .D(LDU.l) .VCLDV.l) ,WORK(l) 


SSVDC  IS  A  SUBROUTINE  TO  REDUCE  A  REAL  NXP  MATRIX  X  BY 
ORTHOGONAL  TRANSFORMATIONS  U  AND  V  TO  DIAGONAL  FORM.  THE 
DIAGONAL  ELEMENTS  S(I)  ARE  THE  SINGULAR  VALDES  OF  X.  THE 
COLUMNS  OF  U  ARE  THE  CORRESPONDING  LEFT  SINGULAR  VECTORS, 

AND  THE  COLUMNS  OF  V  THE  RIGHT  SINGULAR  VECTORS. 

ON  ENTRY 

X  REAL ( LDX, P) .  WHERE  LDX.GE.N. 

X  CONTAINS  THE  MATRIX  WHOSE  SINGULAR  VALUE 
DECOMPOSITION  IS  TO  BE  COMPUTED.  X  IS 
DESTROYED  BY  SSVDC. 

LDX  INTEGER. 

LDX  IS  THE  LEADING  DIMENSION  OF  THE  ARRAY  X. 

N  INTEGER. 

N  IS  THE  NUMBER  OF  COLUMNS  OF  THE  MATRIX  X. 

P  INTEGER. 

P  IS  THE  NUMBER  OF  ROWS  OF  THE  MATRIX  X. 

LDD  INTEGER. 

LDD  IS  THE  LEADING  DIMENSION  OF  THE  ARRAY  0. 
(SEE  BELOW) . 

LDV  INTEGER. 

LDV  IS  THE  LEADING  DIMENSION  OF  THE  ARRAY  V. 
(SEE  BELOW) . 

WORK  REAL(N). 

WORK  IS  A  SCRATCH  ARRAY. 

JOB  INTEGER. 

JOB  CONTROLS  THE  COMPUTATION  OF  THE  SINGULAR 
VECTORS.  IT  HAS  THE  DECIMAL  EXPANSION  AB 
WITH  THE  FOLLOWING  MEANING 

A.EQ.O  DO  NOT  COMPUTE  THE  LEFT  SINGULAR 
VECTORS. 
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C  A.EQ.l  RETCSN  IHE  N  LEFT  SINGULA!  VfcCTDES 

C  IN  C. 

C  A.0E.2  IBTCk  THE  FIRST  MINN.P'  SINGIXaR 

C  VECTORS  IN  D. 

C  B.EQ.O  DO  NOT  OOMPVTE  THE  RIGHT  MNGLTaJI 

C  VECTORS. 

C  B.EQ.l  RETURN  THE  RIGHT  SINGULAR  VE('^>I^ 

C  IN  V. 

C 

C  (m  RETURN 
C 

C  S  REALdOl),  fHERE  lOHMIN  (N*^  I ,  P )  . 

C  THE  FIRST  MININ.?)  ENTRIES  OF  S  CONTAIN  IHE 

C  SINGULAR  VALUES  OF  I  ARRANGED  IN  DESCENDING 

C  ORDER  OF  MAGNITUDE. 

C 

C  B  REALIP). 

C  B  0RDINARIL7  CONTAINS  ZEROS.  BOHEVER  SEE  THE 

C  DISCUSSION  OF  INFO  FOR  EXCEPTIONS. 

C 

C  U  REALILDU.K).  WHERE  LDU.GE.N.  IF  JOBA.EQ.I  THEN 

C  X.EQ.N.  IF  J0BA.GE.2  THEN 

C  K.EQ. MININ, P). 

C  U  CONTAINS  THE  MAIXIl  OF  RIGHT  SINGULAR  VECTORS. 

C  U  IS  NOT  REFERENCED  IF  JOBA.EQ.O.  IF  N.LE.P 

C  OR  IF  J(»A.EQ.2.  THEN  U  MAY  BE  IDENTIFIED  WITH  X 

C  IN  THE  SUBROUTINE  CALL. 

C 

C  V  REALILDV.P).  WHERE  LDV.GE.P. 

C  V  CONTAINS  THE  MATRIX  OP  RIGHT  SINGULAR  VECTORS. 

C  V  IS  NOT  REFERENCED  IF  JGB. m.O.  IF  P.LE.N, 

C  THEN  V  MAI  BE  IDENTIFIED  WITH  X  IN  THE 

C  SUBROUTINE  CALL. 

C 

C  INFO  INTEGER. 

C  THE  SINGULAR  VALUES  (AND  THEIR  CORRESPONDING 

C  SINGULAR  VECTORS)  S(INF(H1),S( INF(H2 S(N) 

C  ARE  CORRECT  (HERE  M-MIN(N,P}).  THUS  IF 

C  INFO.EQ.O.  ALL  THE  SINGULAR  VALUES  AND  THEIR 

C  VECTORS  ARE  CORRECT.  IN  ANY  EVENT,  THE  MATRIX 

C  B  -  IRANS(U)*X«V  IS  THE  BIDIAGONAL  MATRIX 

C  WITH  TEE  ELEMENTS  OF  S  ON  ITS  DIAGONAL  AND  THE 

C  ELEMENTS  OF  E  ON  ITS  SUPER-DIAGONAL  (TRANS(U) 

C  IS  THE  TRANSPOSE  OF  U) .  THUS  THE  SINGULAR 

C  VALUES  OF  X  AND  B  ARE  THE  SANE. 

C 

C  LINPACK.  THIS  VERSION  DATED  03/19/79  . 

C  G.W.  STEWART,  UNIVERSITY  OF  MARYLAND,  ARGONNE  NATIONAL  LAB. 

C 

C  «•••«  USES  THE  FOLLOWING  FUNCTIONS  AND  SUBPROGRAMS. 

C 

C  EXTERNAL  SHOT 

C  BLAS  SAXPY, SDOT, SSCAL. SSWAP, SNRM2 , SROTG 


oonrj  oon  ono 
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C  FORISAN  ABS.AJMAXl.NAXO.MINO, HOD,  SORT 

C 

C  INTERNAL  VARIABLES 

C 

INTEGER  I,ITER.J.J0BD,K.KASE.EK,L,LL,LLS,LM1,LP1,LS.LU,M.MA1IT, 

•  MM,MHl,MPl,NCT,NCrPl,NCD,NRT,NBTPl 
REAL  SDOT.T.R 

REAL  B, C.CS, EL. EMM1.F.G.SNRM2. SCALE. SHIFT, SL.SM.SN.SMMl.Tl, TEST. 

•  ZTEST 
LOGICAL  WANTD.WANTV 

C 

SET  THE  MAXIMDM  NUMBER  OF  ITERATIONS. 

MAIIT  =30 

DETERMINE  WHAT  IS  TO  BE  COMPUTED. 

WANTU  =  .FALSE. 

WANTV  =  .FALSE. 

JOBD  =  MOD(JOB.100)/10 
NCD  =  N 

IF  (JOBD  .GT.  1)  NCD  =  MINO(N.P) 

IF  (JOBD  .NE.  0)  WANTD  =  .TRUE. 

IF  (MOD(JOB,10)  .NE.  0)  WANTV  =  .TRUE. 

REDUCE  X  TO  BIDIAGONAL  FORM,  STORING  THE  DIAGONAL  ELEMENTS 
IN  S  AND  THE  SUPER-DIAGONAL  ELEMENTS  IN  E. 

INFO  =  0 

NCT  =  MIN0(N-1,P) 

NRT  =  MAX0(0,MIN0(P-2.N)) 

LD  =  MAXO(NCT,NRT) 

IF  (LD  .LT.  1)  GO  TO  170 
DO  160  L  =  1,  LD 
LPl  =  L  +  1 

IF  (L  .GT.  NCT)  GO  TO  20 
C 

C  COMPUTE  THE  TRANSFORMATION  FOR  THE  L-TH  COLDMN  AND 

C  PLACE  THE  L-TH  DIAGONAL  IN  S(L) . 

C 

S(L)  =  SNRM2(N-L+1.X(L,L) ,1) 

IF  (S(L)  .EQ.  O.OEO)  GO  TO  10 

IF  (X(L,L)  .NE.  O.OEO)  S(L)  =  SIGN( S(L) .X(L.L) ) 

CALL  SSCAL(N-L+1.1.0E0/S(L),X(L.L),1) 

X(L.L)  =  l.OEO  +  X(L,L) 

10  CONTINUE 

S(L)  =  -S(L) 

20  CONTINUE 

IF  (P  .LT.  LPl)  GO  TO  50 
DO  40  J  =  LPl,  P 

IF  (L  .GT.  NCT)  GO  TO  30 
IF  (S(L)  .EQ.  O.OEO)  GO  TO  30 


\ 
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APPLY  ISE  1SANSF0SMATI(H4. 

T  -  -SD0T(N-L+1,X(L.L),1.1(L,J),1)/X(L,L) 

CALL  SAXPT(N-L->-l.T.X(L.L),l.X(L.J).l) 

30  CONTINUE 

place  IHB  L-TH  BOV  of  X  INTO  E  FOB  ISE 
SUBSEQUENT  CALCULATION  OF  TBE  BOV  TSANSFOBMATION. 

E(J)  -  X(L,J) 

40  CONTINUE 

50  CONTINUE 

IF  (.NOT.VANIU  .OB.  L  .OT.  NCT)  00  TO  70 

PLACE  THE  TSANSFOBMATION  IN  U  FOB  SUBSEQUENT  BACX 
MULTIPLICATION. 

DO  50  I  -  L.  N 
Ud.L)  -  X(I,L) 

60  CONTINUE 

70  CONTINUE 

IF  (L  .6T.  NBT)  GO  TO  150 

COMPOTE  THE  L-TH  BOV  TSANSFOBMATION  AND  PLACE  THE 
L-TH  SUPBB-DIAOONAL  IN  S(L) . 

E(L)  -  SNBM2(P-L.E(LPl).l) 

IF  (E(L)  .EQ.  O.OEO)  GO  TO  80 

IF  (E(LPl)  .NE.  O.OEO)  E(L)  -  SION(E(L) .E(LPl) ) 

CALL  SSCAL(P-L.1.0EO/E(L},E(LP1) ,1) 

E(LPl)  -  l.OEO  +  E(LPl) 

80  CONTINUE 

E<L)  -  -E(L) 

IF  (LPl  .6T.  N  .OS.  E(L)  .EQ.  O.OEO)  GO  TO  120 

APPLY  THE  TSANSFOBMATION. 

DO  90  I  -  LPl,  N 
VOBK(I)  -  O.OEO 
90  CONTINUE 

DO  100  J  -  LPl.  P 

CALL  SAXPY(N-L.E(J).X(LP1,J),1.V0BX(LP1),1) 

100  CONTINUE 

DO  110  J  -  LPl.  P 

CALL  SAXPY(N-L.-EU)/E(LP1)  .VOBK(LPl)  ,1,X(LP1.J)  .1) 
110  CONTINUE 

120  CONTINUE 

IF  (.NOT.VANTV)  GO  TO  140 

PLACE  TBE  TBANSFOBMATUm  IN  Y  FOB  SUBSEQUENT 
BACX  MULTIPLICATION. 


non  non 
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DO  130  I  -  LPl,  P 
V(I.L)  -  Ed) 

130  CONTINDE 

140  CONTINUE 

150  CONTINUE 
160  CONTINUE 
170  CONTINUE 

SET  UP  IHE  FINAL  BIDIAGONAL  NA7SIX  OR  ORDER  N. 

H  =  MIN0(P,N+1) 

NCTPl  *  NCT  +  1 
NRTPl  =«  NRT  +  1 

IF  (NCT  .LT.  P)  S(NCTPl)  =  X( NCTPl, NCTPl) 

IF  (N  .LT.  M)  S(M)  *  O.OEO 

IF  (NRTPl  .LT.  M)  E(NRTPl)  «  X(NRTP1.»() 

E(M)  =  O.OEO 

IF  REQUIRED.  GENERATE  U. 

IF  (.NOT.WANTU)  GO  TO  300 

IF  (NCU  .LT,  NCTPl)  GO  TO  200 
DO  190  J  =  NCTPl.  NCU 
DO  180  I  -  1.  N 
Ud.J)  =  O.OEO 
180  CONTINUE 

Ud.J)  «  l.OEO 
190  CONTINUE 

200  CONTINUE 

IF  (NCT  .LT.  1)  GO  TO  290 
DO  280  LL  =  1.  NCT 
L  =  NCT  -  LL  +  1 
IF  (S(L)  .EQ.  O.OEO)  GO  TO  250 
LPl  =  L  +  1 

IF  (NCU  .LT.  LPl)  GO  TO  220 
DO  210  J  =  LPl.  NCU 

T  -  -SDOT(N-L+l.U(L,L).l,U(L.J),l)/D(L.L) 
CALL  SAXPY(N-L+l.T.O(L.L).l.D(L.J).l) 

210  CONTINDE 

220  CONTINDE 

CALL  SSCAL(N-L+l.-1.0E0.U(L.L).l) 

U(L,L)  -  l.OEO  +  D(L.L) 

LMl  *  L  -  1 

IF  (LMl  .LT.  1)  GO  TO  240 
DO  230  I  >  1.  LNl 
D(I.L)  »  O.OEO 
230  CONTINUE 

240  CONTINUE 

GO  TO  270 
250  CONTINDE 

DO  260  I  «  1.  N 
D(I.L)  =  O.OEO 
260  CONTINUE 


ooooonoo  n  nnnnnnnn  n  n  e»  non 
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U(L.L)  -  l.OBO 
270  OONTIMIIB 

280  amriNDB 

290  OONTINIIB 

300  OONTINDE 

IF  IT  IS  SBQUIKED.  GENEXAIE  V. 

IF  (.NOT.WANTV)  00  TO  350 

DO  340  LL  «  1.  P 
L  -  P  -  LL  ♦  1 
LPl  -  L  +  1 

IF  (L  .OT.  NBT)  GO  TO  320 
IF  (E(L)  .BO.  O.OBO)  GO  TO  320 
DO  310  J  -  LPl,  P 

T  -  -SD0T(P-L.V(LP1,L),1.V(LP1,J),1)/V(LP1,L) 

CALL  SAXPI(P-L,T,V(LP1,L).1,V(LP1,J).1) 

310  CONTINDB 

320  OONTINDE 

DO  330  I  -  1,  P 
V(I.L)  -  O.OEO 
330  OONTZNOE 

V(L.L)  -  l.OEO 
340  OONTINOE 
330  OONTINDE 

MAIN  ITBBATION  LOOP  FOE  THE  S1N6DLAE  VALDES. 

MM  «  M 

I7EE  -  0 
360  CWTINDE 

QDIT  IF  ALL  IBE  SIN6DLAE  VALDES  HAVE  BEEN  FODND. 

...EXIT 

IF  (M  .BO,  0)  60  TO  620 

IF  TOO  MANY  IlEBATIONS  HAVE  BEEN  FEEPOBMED.  SET 

FLAG  AMD  BEIUSN. 

IF  (ITEM  .LT.  MAXIT)  GO  TO  370 
INFO  -  M 

. EXIT 

GO  TO  620 
370  OONTINDE 

THIS  SECTION  OF  THE  FEOGEAM  INSPECTS  POE 

NEGLIGIBLE  ELEMENTS  IN  THE  S  AND  B  AEEATS.  ON 

OOMFLETKW  THE  VARIABLES  EASE  AND  L  ABE  SET  AS  FCHJLOWS. 

EASE  -  1  IF  S(M)  AND  E(l^l)  ABE  NEGLIGIBLE  AND  L.LT.M 
EASE  -2  IF  S(L)  IS  NEGLIGIBLE  AND  L.LT.M 
EASE  -  3  IF  E(L-l)  IS  NEGU6IBLE,  L.LT.M.  AND 
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S(L).  ....  S(M)  ARE  NOT  NEGLIGIBLE  (OR  STEP). 
EASE  *  4  IF  E(M-l)  IS  NEGLIGIBLE  (CONVERGENCE). 

DO  390  LL  =  1,  N 
L  =  M  -  LL 
...EXIT 

IF  (L  .EQ.  0)  GO  TO  400 
TEST  *  ABS(S(L))  +  ABS(S(L+1)) 

ZTEST  =  TEST  +  ABS(E(L)) 

IF  (ZTEST  .NE.  TEST)  GO  TO  380 
E(L)  »  O.OEO 

. EXIT 

GO  TO  400 


380 

CONTINDE 

390 

CONTINUE 

400 

CONTINDE 

IF  (L  .NE.  M  -  1)  GO 
EASE  =  4 

GO  TO  480 

TO  410 

410 

CONTINUE 

LPl  =  L  +  1 

MPl  =■  M  +  1 

DO  430  LLS  =  LPl, 

HPl 

LS  =  M  -  LLS  + 
...EXIT 

LPl 

IF  (LS  .EQ.  L) 
TEST  =  O.OEO 

GO  TO  440 

IF  (LS  .NE.  N) 

TEST  =  TEST  +  ABS(E(LS)) 

IF  (LS  .NE.  L  +  1)  TEST  =  TEST  +  ABS(E(LS-1) 

ZTEST  =  TEST  + 

ABS(S(LS)) 

IF  (ZTEST  .NE. 

TEST)  GO  TO  420 

S(LS)  =  O.OEO 

. EXIT 

GO  TO  440 

420 

CONTINUE 

430 

CONTINDE 

440 

CONTINUE 

IF  (LS  .NE.  L)  GO 
EASE  =  3 

GO  TO  470 

TO  450 

450 

CONTINUE 

IF  (LS  .NE.  M)  GO 
EASE  «  1 

GO  TO  470 

TO  460 

460 

CONTINDE 

EASE  =  2 

L  =  LS 

470 

CONTINUE 

480 

CONTINDE 

L  *  L  +  1 

PERFORM  THE  TASE  INDICATED  BT  EASE. 

lOf 


GO  10  (490,520.540*570),  CASE 
C 

C  DEFLATE  NEGLIGIBLE  S(M) . 

C 

490  0(»4T1NUB 

MMl  -  M  -  1 
P  -  E(lf-l) 

B(ll-l)  «  O.OEO 
DO  510  EX  -  L.  MKl 
E  -  Mm  -  EE  L 
n  -  S(E) 

CALL  SBOXG(Tl.F.CS.SN) 

S<E)  -  n 

IF  (X  ,BQ.  L)  GO  TO  500 
F  -  -SN«B(E-1) 

B(E-l)  -  CS*E(E-1) 

500  OONTINOE 

IF  (lAMTV)  CALL  SS0T(P.V(1.E) .1,V(1 ,M) .l.CS, SN) 

510  OONTINDE 

GO  TO  610 
C 

C  SPLIT  AT  NBSLIGIBLB  S(L} . 

C 

520  OONTINOE 

F  -  B(W) 

B(W)  -  O.OEO 
DO  530  E  •  L,  M 
n  -  S(E) 

CALL  SBOTGdl.F.CS.SN) 

S(E)  >  n 
F  -  -SN*E(E) 

E(E)  -  CS*E(E) 

IF  (WANTU)  CALL  SS0T(N,0(1.E) ,1.D(1,L-1) .l.CS.SN) 

530  OONTINDE 

GO  TO  610 
C 

C  PBBFOBM  ONE  OB  STEP. 

C 

540  OONTINDE 

C  _ 

C  CALCOLAIB  THE  SHIFT. 

C 

SCALE  -  AMAX1(ABS(S(M)).AB8(S(1»>1)),ABS(E(1(-1)).ABS(S(L)), 
•  ABS(E(L))) 

SM  -  S(M) /SCALE 
SMm  -  S()t-1) /SCALE 
EMm  •  E(M>1) /SCALE 
SL  -  S(L)/SCALE 

Bt  m  R/T  1  /RTAt  R 

B  -  ((SMm  SM)*(8Mm  -  SM)  ENm**2)/2.0E0 
C  -  (SM*EMm)*«2 
SBIFT  "*  0,0B0 

IF  (B  .so!  O.OEO  .AND.  C  .BO.  O.OEO)  GO  TO  550 


non  ononrsrs  non 
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SHIFT  =  SQRT(B**2+C) 

IF  (B  .LT.  O.OEO)  SHIFT  =  -SHIFT 
SHIFT  =  C/(B  +  SHIFT) 

550  CONTINUE 

F  =  (SL  +  SM)*(SL  -  SM)  -  SHIFT 
G  =  SL*EL 

CHASE  ZEROS. 

MMl  °  M  -  1 
DO  560  K  ===  L.  MMl 

CALL  SROTG(F.G.CS.SN) 

IF  (K  .NE.  L)  E(K-l)  *  F 
F  =  CS*S(K)  +  SN*E(K) 

E(K)  =  CS*E(K)  -  SN*S(K) 

G  =  SN*S(K+1) 

S(K+1)  =  CS*S(K+1) 

IF  (WANTV)  CALL  SROT(P.V(l,K) .1 ,V(1 ,K+1) .1 .CS, SN) 
CALL  SROTG(F.G.CS.SN) 

S(K)  =  F 

F  =  CS*E(K)  +  SN*S(K+1) 

S(K+1)  =  -SN*E(K)  +  CS*S(K+1) 

G  =  SN*E(K+1) 

E(K+1)  =  CS*E(K+1) 

IF  (WANTD  .AND.  K  .LT.  N) 

•  CALL  SR0T(N,D(1,K),1,D(1.K+1),1,CS,SN) 

560  CONTINUE 

E(M-l)  =  F 
ITER  =  ITER  +  I 
GO  TO  610 

CONVERGENCE. 

570  CONTINUE 

MAKE  THE  SINGULAR  VALUE  POSITIVE. 

IF  (S(L)  .GE.  O.OEO)  GO  TO  580 
S(L)  *  -S(L) 

IF  (WANTV)  CALL  SSCAL(P,-1.0E0.V(1,L) .1) 

580  CONTINUE 

ORDER  THE  SINGULAR  VALUE. 

590  IF  (L  .EQ.  MM)  GO  TO  600 

C  ...EXIT 

IF  (S(L)  .GE.  S(L+1))  GO  TO  600 
T  -  S(L) 

S(L)  -  S(L+1) 

S(L+1)  -  T 

IF  (WANTV  .AND.  L  .LT.  P) 

•  CALL  SSWAP(P,V(l,L),l,V(l,L+l).l) 

IF  (WANTU  .AND.  L  .LT.  N) 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


CALL  SSfAP(N.D(l,L),l.n(l.L+l).l) 
L  -  L  f  1 


GO  TO  590 

600 

CONTINUE 

ITER  «  0 

M  -  M  -  1 

610 

CONTINUE 

60  TO  3<0 
620  CONTINUE 
BE1TJBN 
END 

SEAL  FUNCTION  SNRH2  (  N.  SX.  INCX) 

INTEOES  NEXT,  N  .  INCX 

INTEOES  NN  .  I  .  J 

BEAL  SX(1).  CUILO.  CUTEI.  HITEST,  SUN,  XHAX.  ZERO,  ONE 
DATA  ZERO.  ONE  /O.OEO.  l.OEO/ 

EUCLU^AN  norm  of  tee  N-VECTOR  stored  in  SXO  with  ST0RA6E 
INCREMENT  INCX  . 

IF  N  .LB.  0  RETURN  WITH  RESULT  -  0. 

IF  N  .GE.  I  THEN  INCX  MUST  BE  .6B.  1 

C.L.LAWSm.  1978  JAN  08 

FOUR  PHASE  METHOD  USING  TWO  BUILT-IN  CONSTANTS  THAT  ARE 
HOPEFULLY  APPLICABLE  TO  ALL  MACHINES. 

CUILO  *•  MAXIMUM  OF  SaRT(U/EPS)  OVER  ALL  KNOWN  MACHINES. 

COTHI  *  MINIMUM  OF  SQItl(V)  OVER  ALL  KNOWN  MACHINES. 

WHERE 

EPS  «  SMALLEST  NO.  SUCH  THAT  EPS  +  1.  .OT.  1. 

U  -  SMALLEST  POSITIVE  NO.  (UNDERFLOW  LIMIT) 

V  -  LARGEST  NO.  (OVERFLOW  LIMIT) 

BRIEF  OUTLINE  OF  ALGORITHM. . 

PHASE  1  SCANS  ZERO  COMPONENTS. 

MOVE  TO  PHASE  2  WHEN  A  COMPONENT  IS  NONZERO  AND  .LE.  CUTLO 

MOVE  TO  PHASE  3  WHEN  A  COMPONENT  IS  .OT.  CUILO 

MOVE  TO  PHASE  4  WHEN  A  COMPONENT  IS  .GE.  CUIHI/M 

WHERE  M  «  N  FOR  XO  REAL  AND  M  »  2*N  FOR  COMPLEX. 

VALDES  FOR  CUTLO  AND  CUTEI.. 

FROM  THE  ENVIRONMENTAL  PARAMETERS  LISTED  IN  THE  IMSL  CONVERTER 
DOCUMENT  THE  LIMITING  VALUES  ARE  AS  FOLLOWS. . 

CUTLO,  S.P.  U/EPS  -  2««(-102)  FOR  HONEYWELL.  CLOSE  SECONDS  ARE 
UNIVAC  AND  DEC  AT  2**(-103) 

THUS  CUTLO  -  2**(-5l)  -  4.44089B-16 
CDTHI,  S.P.  V  -  2*«127  FOR  UNIVAC,  HONEYWELL.  AND  DEC. 

THUS  CUTHI  «  2««(63.5)  -  1.30438B19 
CUTLO,  D.P.  U/EPS  -  2««(-67)  FDR  HONEYWELL  AND  DEC. 

THUS  CUTLO  -  2*«(-33.5)  -  8.23181D-11 
CDTHI.  D.P.  SAME  AS  S.P.  CDTHI  -  1.30438D19 
DATA  CUTLO,  CDTHI  /  8.232D-11.  1.304D19  / 


li 


X 
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C  DATA  CDTLO.  CDTHI  /  4.441E-16.  1.304E19  / 

DATA  CDTLO.  CDTHI  /  4.441E-16.  1.304E19  / 

C 

IF(N  .GT.  0)  GO  TO  10 
SNRM2  °  ZERO 
GO  TO  300 
C 

10  ASSIGN  30  TO  NEXT 
SDN  °  ZERO 
NN  =  N  •  INCX 

C  BEGIN  MAIN  LOOP 

1  =  1 

20  GO  TO  NEXT, (30,  50.  70.  110) 

30  IF(  ABSCSXd))  .GT.  CDTLO)  GO  TO  85 
ASSIGN  50  TO  NEXT 
XMAX  =  ZERO 

PHASE  1.  SDN  IS  ZERO 

50  IF(  SX(I)  .EQ.  ZERO)  GO  TO  200 

IF(  ABS(SX(I))  .GT.  CDTLO)  GO  TO  85 

PREPARE  FOR  PHASE  2. 

ASSIGN  70  TO  NEXT 
GO  TO  105 

PREPARE  FOR  PHASE  4. 

100  I  =  J 

ASSIGN  110  TO  NEXT 
SDN  =  (SDN  /  SX(I))  /  SX(I) 

105  XMAX  =  ABS(SX(I)) 

GO  TO  115 

PHASE  2.  SDN  IS  SMALL. 

SCALE  TO  AVOID  DESTRDCTIVE  UNDERFLOW. 

70  IF(  ABS(SX(I))  .GT.  CDTLO  )  GO  TO  75 

COMMON  CODE  FOR  PHASES  2  AND  4. 

IN  PHASE  4  SUM  IS  LARGE.  SCALE  TO  AVOID  OVERFLOW. 

110  IF(  ABS(SX(I))  .LE.  XMAX  )  GO  TO  115 

SDN  =  ONE  SDM.«  (XMAX  /  SX(I))**2 
XMAX  =  ABS(SX(I)) 

00  TO  200 

115  SDM  =  SDM  +  (SX(I)/XMAX)**2 
GO  TO  200 


PREPARE  FOR  PHASE  3. 
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75  SDH  -  (SDH  •  IMAX)  •  ZMAX 


FOE  SEAL  OR  D.P.  SET  HITEST  -  CnTSl/N 
FOB.  OQliFLEX  SET  HITEST  -  CinSI/(2*N) 

85  HITEST  -  CDTHl/FLOAK  N  ) 

PHASE  3.  SDM  IS  MID-RANGE.  NO  SCALING. 

DO  95  I  -I.NN.INCX 
IF(ABS(SX(J))  .GE.  HITEST)  GO  TO  100 
95  SDM  -  SOM  SX(J}**2 
SNRM2  «  9(UET(  SDM  ) 

GO  TO  300 

200  CONTINDB 

I  -  I  ♦  INCX 

IF  (  I  .LE.  NN  )  GO  TO  20 
END  OF  MAIN  LOOP. 

COMPOTE  SODARE  ROOT  AND  ADJDST  FOR  SCALING. 

SNRM2  «  2MAZ  *  SQRT(SOM) 

300  CONTINOB 
RETURN 
END 

SDBROOTINE  SSCALCN.  SA.  SX.INCX) 

SCALES  A  VECTOR  HI  A  CONSTANT. 

USES  UNROLLED  LOOPS  FOR  INCREMENT  EQUAL  TO  1. 

JACK  DCmOARRA.  LINPACX.  3/11/78. 

REAL  SA.SZ(l) 

INTEGER  I.INCX.M.MPl.N.NINCX 

IF(N.LE.O)RETURN 
IF(INCX.BQ.l)GO  TO  20 

com  FOR  INCREMENT  NOT  EQUAL  TO  1 

NINCX  -  N*IMCX 
DO  10  I  -  I.NINCX.INCX 
SXd)  -  SA«SX(I) 

10  CONTINUE 
RETURN 

com  FOR  INCREMENT  EQUAL  TO  1 


J 


CLEAN-UP  LOW 
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20  M  =  M0D(N,S) 

1F(  M  .EQ.  0  }  GO  TO  40 
DO  30  I  =  l.M 
SX(I)  =  SA*SX(I) 

30  CONTINUE 

IF(  N  .LT.  5  )  REIU8N 
40  MPl  =  M  +  1 

DO  50  I  =  MP1.N.5 
SX(I)  =  SA*SX(I) 

SX(I  +  1)  =  SA*SX(I  +  1) 

SX(I  +  2)  =  SA*SX(I  +  2) 

SX(I  +  3)  =  SA*SX(I  +  3) 

SX(I  +  4)  =  SA*SX(I  +  4) 

50  CONTINUE 
RETURN 
END 

REAL  FUNCTION  SD0T(N. SX. INCX, SY. INCY) 

FORMS  THE  DOT  PRODUCT  OF  TWO  VECTORS. 

USES  UNROLLED  LOOPS  FOR  INCREMENTS  EQUAL  TO  ONE. 
JACK  DONGARRA.  LINPACK,  3/11/78. 

REAL  SX(1) ,SY(1) .STEMP 
INTEGER  I.INCX.INCY.IX.IY.M.MPl.N 

STEMP  =  O.OEO 
SDOT  =  O.OEO 
IF(N.LE.O)RETURN 

IFdNCX. EQ.l. AND. INCY. EQ.DGO  TO  20 

CODE  FOR  UNEQUAL  INCREMENTS  OR  EQUAL  INCREMENTS 
NOT  EQUAL  TO  1 

IX  =  1 
lY  =  1 

IF(INCX.LT.O)IX  =  (-N+1)*INCX  +  1 
IF(INCY.LT.O)IY  =  (-N+1)*INCY  +  1 
DO  10  I  =  l.N 

STEMP  =  STEMP  +  SX( IX) •SYC lY) 

IX  =  IX  +  INCX 

lY  =  lY  +  INCY 

10  CONTINUE 
SDOT  =  STEMP 
RETURN 
C 

C  CODE  FOR  BOTH  INCREMENTS  EQUAL  TO  1 

C 

C 

C  CLEAN-UP  LOOP 

C 

20  M  =  MOD(N.5) 

IF(  M  .EQ.  0  )  GO  TO  40 
DO  30  I  -  l.M 
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SIEMP  =  STEMP  +  SX(I)*SY(I) 

30  CONTINUE 

IF(  N  .LT.  5  )  GO  TO  60 
40  MPl  =  M  +  1 

DO  50  I  =  MP1.N.5 

STEMP  =  STEMP  +  SX{I)*SY(I)  +  SX(I  +  1)*SY(I  +  1)  ♦ 

•  SKI  +  2)*SY(I  +  2)  +  SKI  +  3)*SY(I  +  3)  +  SKI  4)*SY(I  +  4) 
50  CONTINUE 
60  SDOT  =  STEMP 
RETURN 
END 

SUBROUTINE  SAXPY(N. SA. SX. INCX. SY. INCY) 

C 

C  CONSTANT  TIMES  A  VECTOR  PLUS  A  VECTOR. 

C  USES  UNROLLED  LOOP  FOR  INCREMENTS  EQUAL  TO  OM 
C  JACK  DONGARRA.  LINPACK,  3/11/78. 

C 

REAL  SX(1) ,SY(1) .SA 
INTEGER  I.INCX.INCY.IX.IY.M.MPl.N 

IF{N.LE.O) RETURN 
IF  (SA  .EQ.  0.0)  RETURN 
IF(INCX.EQ.l.AND.INCY.EQ.l)GO  TO  20 

CODE  FOR  UNEQUAL  INCREMENTS  OR  EQUAL  INCREMENTb 
NOT  EQUAL  TO  1 

IX  =  1 
lY  =  1 

IF(INCX.LT.O)IX  =  (-N+1)*INCX  +  1 
IF(INCY.LT.O)IY  =  (-N+1)*INCY  +  1 
DO  10  I  =  1,N 

SY(IY)  =  SY(IY)  +  SA*SX(IX) 

IX  =  IX  +  INCX 

lY  =  lY  +  INCY 

10  CONTINUE 
RETURN 

C  CODE  FOR  BOTH  INCREMENTS  EQUAL  TO  1 

C 

C 

C  CLEAN-UP  LOOP 

C 

20  M  >  MOD(N.4) 

IF(  M  .EQ.  0  )  GO  TO  40 
DO  30  I  «  l.M 

SY(I)  =  SYd)  +  SA*SX(I) 

30  CONTINUE 

IF(  N  .LT.  4  )  RETURN 
40  MPl  -  M  +  1 

DO  50  I  =  MP1.N.4 

ST(I)  =  SYd)  +  SA*SXd) 

SYd  +  1)  »  SYd  +  1)  +  SA'SXCI  +  1) 
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SY(I  +2)  =  SY(I  +  2)  +  SA*SX(I  +  2) 

SY(I  +  3)  =  SY(I  +  3)  +  SA*SX(I  +  3) 

50  CONTINUE 
RETURN 
END 

SUBROUTINE  SR0TG( SA. SB.C. S) 

CONSTRUCT  GIVENS  PLANE  ROTATION. 

JACK  DONGARRA.  LINPACI.  3/11/78. 

REAL  SA. SB, C.S, ROE, SCALE, R.Z 

ROE  =  SB 

IF(  ABS(SA)  .GT.  ABS(SB)  )  ROE  =  SA 
SCALE  =  ABS(SA)  ABS(SB) 

IF(  SCALE  .NE.  0.0  )  GO  TO  10 
C  =  1.0 
S  =  0.0 
R  =  0.0 
GO  TO  20 

10  R  =  SCALE*SQRT((SA/ SCALE) **2  +  { SB/SCALE) **2) 

R  =  SIGN(1.0.ROE)*R 
C  =  SA/R 
S  =  SB/R 
20  Z  =  1.0 

IF(  ABS(SA)  .GT.  ABS(SB)  )  Z  =  S 

IF(  ABS(SB)  .GE.  ABS(SA)  .AND.  C  .NE.  0.0  )  Z  =  1.0/C 

SA  =  R 

SB  =  Z 

RETURN 

END 

SUBROUTINE  SROT  (N. SX. INCX, SY. INCY, C, S) 

APPLIES  A  PLANE  ROTATION. 

JACK  DONGARRA.  LINPACK,  3/11/78. 

REAL  SX(1) .SY(1) .STEMP.C.S 
INTEGER  I.INCX.INCY.IX.IY.N 

IF (N.LE.O) RETURN 

IFdNCX. EQ.l. AND.  INCY. EQ.DGO  TO  20 

CODE  FOR  UNEQUAL  INCREMENTS  OR  EQUAL  INCREMENTS  NOT  EQUAL 
TO  1 

IX  =  1 
lY  =  1 

IF(INCX.LT.O)IX  =  (-N+1)*INCX  +  1 
IF(INCY.LT.O)IY  =  (-N+1)*INCY  +  1 
DO  10  I  =  1,N 

STEMP  =  C*SX(IX)  +  S*SY(IY) 

SY(IY)  =  C*SY(IY)  -  S*SX(IX) 

SX(IX)  =  STEMP 
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II  =  II  +  INC! 
lY  =•  lY  +  INCY 
10  CONTINUE 
RETURN 

CODE  FOR  BOTH  INCREMENTS  EQUAL  TO  1 

20  DO  30  I  =  l.N 

STEMP  =  C*SX(I)  +  S*SY<I) 

SY(I)  »  C*SY(I)  -  S*SX(I) 

SKI)  =■  STEMP 
30  CONTINUE 
RETURN 
END 

SUBROUTINE  SSWAP  (N. SX. INCX, SY. INCY) 

INTERCHANGES  TWO  VECTORS. 

USES  UNROLLED  LOOPS  FOR  INCREMENTS  EQUAL  TO  1. 

JACK  DONGARRA,  LINPACK.  3/11/78. 

REAL  SKI)  .SY(1)  .STEMP 
INTEGER  I.INCX.INCY.IX.IY.M.MPl.N 

IF(N.LE.O)RETURN 

IF(INCX. EQ.l. AND. INCY. EQ.DGO  TO  20 

CODE  FOR  UNEQUAL  INCREMENTS  OR  EQUAL  INCREMENTS  NOT  EQUAL 
TO  1 

IX  =  1 
lY  =  1 

IF(INCX.LT.O)IX  =  (-N+1)*INCX  +  1 
IF(INCY.LT.O)IY  =  (-N+1)*INCY  +  1 
DO  10  I  =  l.N 
STEMP  -  SKIX) 

SKIX)  =«  SY(IY) 

SY(IY)  =«  STEMP 

IX  =  IX  +  INCX 
lY  =  lY  +  INCY 
10  CONTINUE 
RETURN 
C 

C  CODE  FOR  BOTH  INCREMENTS  EQUAL  TO  1 

C 

C 

C  CLEAN-UP  LOOP 

C 

20  M  >  MOD(N.3) 

IF(  M  .BQ.  0  )  GO  TO  40 
DO  30  I  *=  l.M 
STEMP  =  SKI) 

SKI)  “  SY(I) 

SY(I)  =  STEMP 


A 


i 


lis 


30  CONTINOE 

IF(  N  .LT.  3  )  RETURN 
40  MPl  =  M  +  1 

DO  50  I  =  MP1,N.3 
STEMP  =  SX(I) 

SX(I)  =  ST(I) 

SKI)  =  STEMP 
STEMP  =  SX(I  +  1) 

SX(I  +  1)  =  SKI  +  1) 
SKI  +  1)  =  STEMP 
STEMP  =  SX(I  +  2) 

SX(I  +  2)  =  SKI  +  2) 
Sr(I  +  2)  =  STEMP 
50  CONTINDE 
RETURN 
END 
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